NASA CR-2391 



A FINITE ELEMENT SOLUTION ALGORITHM 
FOR THE NAVIER-STOKES EQUATIONS 


by A. J. Baker 


Prepared by 

BELL AEROSPACE COMPANY 

Buffalo, N.Y. 14240 

Jar Langley Research Center 





NATIONAL AERONAUTICS AND SPACE ADMINISTRATION • WASHINGTON, D. C. 


JUNE 1974 





1. Report Ho. 2. Government Accession No. 

NASA CR-2391 

4. Title and Subtitle 

A Finite Element Solution Algorithm for the 
Navier-Stokes Equations 

7. Author(s) 

A.J. Baker 

9. Performing Organization Name and Address 

Bell Aerospace Company 
P.O. Box One 
Buffalo, New York 14240 

12. Sponsoring Agency Nome and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 


3- Recipient’s Catalog No. 

5. Report Date 

June 1974 

6. Performing Organization Code 

8. Performing Organization Report No. 

D9 198-9 50001 

10. Work Unit Ne. 

11. Contract or Grant No. 

NAS 1-11809 

13. Type of Report and Period Covered 

Contractor Report 

14. Sponsoring Agency Code 


15. Supplementary Notes 

Final Report 

16. Abstract 

A Finite element solution algorithm is established for the two-dimensional Navier-Stokes equa- 
tions governing the steady-state kinematics and thermodynamics of a variable viscosity, compressible 
multiple-species fluid. For an incompressible fluid, the motion may be transient as well. The primi- 
tive dependent variables are replaced by a vorticity-streamfunction description valid in domains 
spanned by rectangular, cylindrical and spherical coordinate systems. Use of derived variables pro- 
vides a uniformly elliptic partial differential equation description for the Navier-Stokes system, and foi 
which the Finite element algorithm is established. Explicit non-linearity is accepted by the theory, 
since no psuedo-variational principles are employed, and there is no requirement for either computa- 
tional mesh or solution domain closure regularity: Boundary condition constraints on the normal flux 
and tangential distribution of all computational variables, as well as velocity, are routinely piecewise 
enforceable on domain closure segments arbitrarily oriented with respect to a global reference frame. 

A natural coordinate function description is established for the Finite element approximation 
functions that renders matrix evaluation straightforward and accurate. The several explicitly non- 
linear terms in the equation system are shown to be expressible in terms of standard matrix forms. 

The consequence of integrating various terms by parts is thoroughly evaluated. The COMOC com- 
puter program system embodies the established finite element algorithm and matrix generator package 
built upon the natural function description. Numerical solutions for incompressible internal flow 
problems have illustrated solution accuracy, convergence, and versatility. Flow Fields with imbedded 
regions of recirculation have been computed as well as transient solution domain closure shapes. 
•Evaluation of pressure computations has been performed as well as extension to compressible flow 
field solutions. 


17. Key Words (V 'ected by Author(s)) 

Finite Element 
Navier-Stokes Equations 
Natural Coordinate Functions 
Computer Program 


18. Distribution Statement 


Unclassified - Unlimited 


19. Security Clossif. (of this report) 

20. Security Classif. (of this page) 

21 . No. of Pages 

22. Price* 

Unclassified 

Unclassified 

75 

$ 3.75 


For sale by the National Technical Information Service, Springfield, Virginia 22151 



Page Intentionally Left Blank 



CONTENTS 


Page 

SUMMARY 1 

INTRODUCTION 1 

NOMENCLATURE • 3 

THE NAVIER-STOKES EQUATIONS FOR A MULTICOMPONENT VISCOUS FLUID . . 6 

DIFFERENTIAL EQUATION DEVELOPMENT FOR THE DESIRED PROBLEM 

CLASS IN FLUID MECHANICS 8 

Continuity Equation 8 

Compatibility Equation 10 

The Curl of the Navier-Stokes Equation • 10 

The x 3 Component of the Navier-Stokes Equation 13 

The Species Continuity Equation 14 

The Energy Equation 14 

Recovery of Pressure 16 

FINITE ELEMENT SOLUTION ALGORITHM FOR THE NAVIER-STOKES 

EQUATIONS 18 

FINITE ELEMENT MATRIX GENERATION 22 

Planar Finite Elements for the Two-Dimensional Navier-Stokes Equations 25 

Planar Ring Finite Elements for the Axisymmetric Incompressible Navier-Stokes 

Equations 31 

COMOC COMPUTER PROGRAM ' 34 

NUMERICAL RESULTS 34 

CONCLUDING REMARKS 54 

APPENDIX A: CARTESIAN TENSORS IN EUCLIDEAN SPACE 59 

APPENDIX B: OBSERVATIONS ON THE FINITE ELEMENT SOLUTION 

ALGORITHM FOR THE NAVIER-STOKES EQUATIONS USING 
LINEAR NATURAL COORDINATE APPROXIMATION 

FUNCTIONALS 62 

APPENDIX C: SYSTEM INVARIANCE UNDER COORDINATE 

TRANSFORMATION 69 

REFERENCES 72 


iii 


ILLUSTRATIONS 


Figure Page 

1 Intrinsic Finite Element Domains for Simplex Approximation Functions 23 

2 Establishment of Vorticity Boundary Condition Statement 31 

3 Discretization of Rectangular Duct into 1 44 Triangular Finite Elements 35 

4 Computed Fully Developed Longitudinal Velocity Distributions, Re = 200 ... 36 

5 Computed Fully Developed Streamfunction and Vorticity for Duct Flow, 

Re = 200 : - 37 

6 Longitudinal Velocity Distributions for Incompressible Duct Flow, Re = 200 .. . 38 

7 Computed Transient Temperature Distribution in Axisymmetric Quiescent 

Duct 40 

8 Discretization Influence on Vorticity and Streamfunction Near a Wall 41 

9 Discretization Influence on Computed Streamfunction 42 

10 Flow Over a Rearward Facing Step 43 

1 ] Discretization of Rearward Facing Step in a Rectangular Duct into 21 1 

Triangular Finite Elements 44 

1 2 COMOC Computed Streamfunction and Vorticity Distribution for Transient 

Flow Over a Rearward Facing Step 45 

1 3 COMOC Computed Steady-State Streamfunction and Vorticity for Flow Over 

a Rearward Step in Irregular Shaped Duct, and Over a Rearward Step with 
Internal obstacle, Re = 200 47 

14 Centerplane Pressure Decay in Steady-State Duct Flow, Re = 200 49 

1 5 Centerplane Pressure Distribution for Steady Flow Over a Rearward Facing 

Step, Re = 200 . '...' 50 

16 Flow into a Compression Corner 52 

1 7 Supersonic Flow Over a Rearward Facing Cone-Step 52 

18 Computed Steady-State Distributions for Compressible Flow Over a 

Rearward Step, Re = 200, = 0.6. . : : ■. 55 

19 Finite Element Discretizations for Supersonic Boundary Flows 56 

20 Computed Steady-State Distributions for Supersonic Boundary Layer 

Flow, M^ = 3.0 57 

B-l Two-Dimensional Finite Element 63 

B-2 Adjacent Finite Elements 65 

C-l Coordinate Systems for the Finite Element Solution Algorithm 70 


TABLES 

; I 

Number Page 

1 Coefficients in Generalized Elliptic Differential Equation Description, 

Equation (76) 20 

2 Implicit Definition of Simplex Natural Coordinate Functions 24 

3 Integrals of Natural Coordinate Function Products Over Finite Element 

Domains 24 

4 Standard Finite Element Matrix Forms for Simplex Functionals in Two- 

Dimensional Space 26 

5 Additional Standard Matrix Forms for Simplex Functionals in Axisymmetric 

Space i 33 


IV 


A FINITE ELEMENT SOLUTION ALGORITHM 
FOR THE NAVIER-STOKES EQUATIONS 
By 

A.J. Baker 

Bell Aerospace Company 
SUMMARY 


A finite element solution algorithm is established for the two-dimensional Navier-Stokes equa- 
tions governing the steady-state kinematics and thermodynamics of a variable viscosity, compressible 
multiple-species fluid. For an incompressible fluid, the motion may be transient as well. The primi- 
tive dependent variables are replaced by a vorticity-streamfunction description valid in domains span- 
ned by rectangular, cylindrical and spherical coordinate systems. Use of derived variables provides a 
uniformly elliptic partial differential equation description for the Navier-Stokes system, and for 
which the finite element algorithm is established. Explicit nonlinearity is accepted by the theory, 
since no psuedo-varjational principles are employed, and there is no requirement for either compu- 
tational mesh or solution domain closure regularity. Boundary condition constraints on the normal 
flux and tangential distribution of all computational variables, as well as velocity, are routinely 
piecewise enforceable on domain closure segments arbitrarily oriented with respect to a global ref- 
erence frame. 

A natural coordinate function description is established for the finite element approximation 
functions that renders matrix evaluation straightforward and accurate. The several explicitly non- 
linear terms in the equation system are shown to be expressible in terms of standard matrix forms. 
The consequence of integrating various terms by parts is thoroughly evaluated. The COMOC com- 
puter program system embodies the established finite element algorithm and matrix generator pack- 
age built upon the natural function description. Numerical solutions for incompressible internal flow 
problems have illustrated solution accuracy, convergence, and versatility. Flow fields with imbedded 
regions of recirculation have been computed as well as transient solution domain closure shapes. 
Evaluation of pressure computations has been performed as well as extension to compressible flow 
field solutions. 


INTRODUCTION 


Numerical solution of a variety of field problems in mechanics has been made possible with 
the advent of the large digital computer. Development of solution procedures for specific disciplines 
has been highly problem oriented. As a result, little cross-fertilization has occurred that takes advan- 
tage of the uniform mathematical description for field problems in continuum mechanics. Specifi- 
cally. finite difference methods have been almost universally employed for computational fluid mech- 
anics, the electric analogy approach has been applied to heat conduction and seepage problems, while 
the finite element method has found wide acceptance for analysis of complex structural systems. Each 
of these approaches can now be viewed in a unifying manner as application of specific criteria within 
the Method of Weighted Residuals (MWR) (Reference 1), and satisfaction in a weighted average sense 
of the governing differential equations written on the discretized equivalent of the pertinent depend- 
ent variables. The select choice of weighting and approximation functions renders each approach identi- 



fiable, including classical integral approaches like von Karman-Pohlhausen and Integral Method of 
Strips in boundary layer flow, and the many variations of Galerkin, Kantorovich and collocation 
methods. 

Over the years, the finite element procedure has proven highly adaptable to solution of linear 
elliptip boundary value problems involving complex boundary conditions applied on irregularly 
shaped, non-coordinate surface solution domain closures. Its adaptation to quasilinear and/or nonlin- 
ear problem classes — for example, solution of the Navier-Stokes equations — has been held in abey- 
ance, at least partially, by the, common misconception that the method was constrained to problems 
having a differential equation description that could be equivalently cast into an extremization prin- 
ciple. This in turn led to derivation of a plethora of psuedo-variational statements for nonstationary 
field problems (see for example Reference 2) in an attempt to render the finite element method direct- 
ly applicable to the broader problem class. Finite element solution of transient linear heat conduction 
in stationary continua drew much attention in the late 1960’s (References 3, 4), and inclusion of 
implicit nonlinearity was eventually acceptable (Reference 5) within the constraints of “local” ex- 
tremization. Extension to the Navier-Stokes equations using these concepts appeared theoretically 
difficult without gross assumptions on the concept of a local potential. However, about 1970 the 
connection between Galerkin criteria within MWR and equivalent extremization principles was 
established, the key items being performance of an integration by parts, identification of a numerical 
Lagrange multiplier, and Boolean assembly of the local finite element matrix equations into the glob- 
al description. It could be readily proved that such a development, for a linear equation, produced a 
computational form that was identical to extremization of the equivalent stationary principle. How- 
ever, since no linearity constraint existed in its derivation, the theory for finite element solution of 
nonlinear equations was established, especially for specific forms of the Navier-Stokes equations 
(References 6-8). 

In this report, a finite element solution algorithm is established for the two-dimensional 
Navier-Stokes equations governing the kinematics and thermodynamics of a variable viscosity, com- 
pressible, multiple-species fluid. The primitive dependent variables are replaced by a vorticity- 
streamfunction description, which provides a uniformly elliptic differential equation system descrip- 
tion for all computational variables. The preferred differential equation systems are established in 
rectangular, cylindrical and spherical Cartesian coordinate systems. The finite element algorithm is 
derived for the generalized, nonlinear elliptic boundary value problem of mathematical physics and 
contains no requirements for either computational mesh or solution domain closure regularity. 
Boundary condition constraints on the normal flux and tangential distribution of all computational 
dependent variables, as well as velocity, are routinely piecewise enforceable on domain closure 
segments arbitrarily oriented with respect to a global reference frame. The intrinsic finite element 
shapes for one-, two-, and three-dimensional domains spanned by linear approximation functions are 
the line, triangle and tetrahedron respectively. The area-coordinate concept of structural mechanics 
has been utilized to establish a natural coordinate function description for the finite element approx- 
imation to the Navier-Stokes equations. Adaptation of these functions represents a significant advance 
that has for the first time allowed analytical formulation of the complexly nonlinear matrix repre- 
sentations for specific terms in the equation system — in particular, the vorticity transport equation. 
The consequence of integration by parts of select nonlinear terms in the differential equation system 
is examined in detail, and the assumptive constraints are established according to which the generated 
surface integrals can be neglected or made to cancel. Establishment of the natural coordinate func- 
tion description has allowed recognition of the many standard matrix forms that constitute. the finite 
element algorithm for the nonlinear terms in the equation system and provides broad insight into the 
mechanics of the algorithm. 
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The finite element solution algorithm for the characteristic equation system has been embod- 
ied into the COMOC (Computational Continuum Mechanics) computer program system, specific var- 
iants of which have produced solutions in three-dimensional subsonic and supersonic viscous flow fields 
(References 9, 10) and nonlinear transient heat conduction (Reference 1 1) in addition to the Navier- 
Stokes solutions to be discussed. The COMOC system consists of four basic macro-modules, the first 
of which is Input wherein discretizations are formed and dependent variables initiated. The Geometry 
module establishes the non-standard element matrices for each finite element of the discretization and 
evaluates the matrix multipliers required for the standard matrices. All computations are performed 
in a local Cartesian coordinate system with automatic accounting for the lack of coordinate trans- 
formation invariance for other than rectangular Cartesian systems. The Integration module embodies 
solution algorithms for large-order ordinary differential and algebraic equation systems, either of 
which is produced by application of the finite element algorithm to the parent partial differential 
equation system. For discretized-equivalent initial value problems, the basic operation is evaluation 
of the derivative vector on an element basis and assembly of the global vector using Boolean algebra. 
For solution of an algebraic system, the diffusion (that is, “stiffness” in elasticity) matrix is formed 
on an element basis and the equation system inhomogeniety evaluated, if present. Each is then 
assembled into the global representation through Boolean algebra. The matrix equation system rank 
is established by evaluating boundary condition constraints, and the effective inverse found by equa- 
tion solver techniques. The final Output module serves its standard function. 

The Navier-Stokes Variant of COMOC has been evaluated for several problems to assess solu- 
tion accuracy, convergence, versatility, and overall performance of the system. Accuracy studies 
were performed for developing flow in a duct, and convergence of velocity with discretization is 
numerically demonstrated. Factors affecting solution accuracy have been evaluated including non- 
uniform discretizations, the vorticity boundary condition, and the use of a “condensed” mass matrix 
for solving transient problems. A basic character of full Navier-Stokes solutions is prediction of im- 
bedded regions of recirculation, and the finite element algorithm is assessed to correctly predict such 
phenomena for a sample problem without resorting to special boundary condition techniques. A 
stable solution for the rearward step problem has been obtained for the largest flow Reynolds number 
reported in the literature. Variations on this problem have illustrated stable solutions for time-depend- 
ent solution domain closure configurations and solutions in multiply-connected domains. Evaluation 
of the accuracy of pressure computations for these problems is presented as well as extension to flow 
fields with variable density. The computed results have yielded a favorable assessment of the viability 
of the solution algorithm and its computational embodiment. 

NOMENCLATURE 


a boundary condition coefficient; unit vector 

b body force unit vector 

B body force 

c expansion coefficient; specific heat 

D diffusion coefficient 

e general alternating tensor 

E symmetric velocity gradient 

Ec Eckert Number 
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function of known argument 
Froude Number 

mass flow; function of known argument 

metric coefficient 

stagnation enthalpy 

rectangular Cartesian unit vectors 

determinant of Cartesian space metric v 

function of known argument; thermal conductivity 

generalized diffusion coefficient 

characteristic length; natural coordinate function; arc length 
Lewis Number 

Mach Number; number of finite elements 
unit normal vector; summation limit 
differential operator 
pressure 
Prandtl Number 

generalized dependent variable; energy source term 

spherical Cartesian coordinate system 

position vector; radius 

domain of space variables ; gas constant 

Reynolds Number 

displacement 

mass source term 

time 

temperature; total stress tensor; tensor 

velocity 

vector 

weighting function 

rectangular Cartesian coordinate system 
independent space variables 
rectangular Cartesian scalar components 
mass fraction 

cylindrical Cartesian coordinate system 

coefficient for spherical coordinates; mass fraction; direction cosine 
ratio of specific heats 
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finite element coefficient matrix 

Kronecker delta 

closure of solution domain 

alternating tensor; strain; vector 

differential operator 

angle 

coefficient 
Lagrange multiplier 
viscosity 

unit normal vector 

density I , 

surface integral kernel 
stress tensor; domain integral kernel 
angle; function 

x 3 scalar component of streamfunction 
vector potential 

x 3 spalar component of vorticity ; wryness coefficient 
vorticity vector 
column matrix 

square matrix 
diagonal matrix 

Superscripts and Subscripts 

approximate solution 
temporal derivative 
matrix transpose 
unit vector 
vector 

average; constrained to domain closure 
tensor indices 

pertaining to m th subdomain (finite element) 
reference condition 



THE NAVIER-STOKES EQUATIONS FOR A MULTICOMPONENT 

VISCOUS FLUID 


The complete description of a fluid dynamical state is contained within the solution of the 
system of coupled, nonlinear, second order partial differential equations describing the local con- 
servation of species, mass, linear (or angular) momentum, and energy, in conjunction with an appro- 
priate specification of constitutive behavior and applicable initial and boundary conditions. In 
Cartesian tensor notation (see Appendix A) the conservation form of these differential equations is 


respectively 

(pY a ), t = -[pUj Y a - pDY a , j] ; j + S a (1) 

P, t = - tP u i ) . j (2) 

(„ u.), , - - 1* u, Uj + pSjj • r # ] . j +pBj . ( 3 ) 

(pH - p) , t = - [pu. H - r..u. - kT, .] . . + p Q (4) 


The dependent variables have their usual fluid dynamic interpretation, with p the mass den- 
sity, Y a the mass fraction of the a** 1 species, uj the local velocity vector, S a the generation term of 
the a th species, Bj the applicable body force, r y the viscous stress tensor, H the stagnation enthalpy, 
p the pressure, T the temperature, and Q the local heat generation rate. In equations (1) - (4) the 
comma denotes partial differentiation of a scalar, while the semicolon indicates the Vector derivative 
of a Cartesian tensor, the Cartesian equivalent of covariant differentiation of the contravarient com- 
ponents of a general tensor. The comma by t indicates the partial derivative with respect to time. 

The solution of the equation system (1) - (4) requires specification of constitutive relation- 
ships between the dependent variables and the diffusion coefficient D, the stress tensor Ty, the 
viscosity p and thermal conductivity k. A specification of an equation of state relating the thermo- 
dynamic variables is also required. For the laminar flow of a Newtonian fluid, the dynamic relations 


are contained within Stokes’ viscosity law 

T ij = 2p[Eij-y E kk Sij] (5) 

Eij=y (ui ;j + uj -i) (6) 

The thermodynamic properties are typically 

p = p(p,Y?T) (7) 

H = H (p, T, ui) = [I p dT (8) 

0 
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D a 

= D(Y a , Uj) 

(9) 

t 

11 

s 

H 

w 

(12) 

S a 

= S(Y a ) 

(10) 

Q - Q(Y a , T) 

(13) 

k 

= k (T) 

(ID 




The functional dependence indicated in equations (7) - (13) can be determined experimen- 
tally. Equations (1) - (13) form a deterministic system for the dependent variables Y a , p, uj, and T, 
and represent a well posed, initial-boundary value problem in mathematical physics for proper speci- 
fication of these constraints. 


Before transforming equations (1) - (13) into the desired form, it is convenient to nondimen- 
sionalize all variables to extract the useful non-dimensional groupings, the magnitudes of which 
characterize various flow regimes. Denoting non-dimensional variables by an asterisk. let 


>< 

II 

r 

>< 

(14) 

T = ToqT* 

(20) 

u i = U^u* 

(15) 

C P = C Poo C P* 

(21) 

II 

(16) 

k = 

Pr 

(22) 

P - PooP* 

(17) 

^ kocTo° 

q = ^ q * 

(23) 

P = PooU 2 ooP * 

M = MooP* 

(18) 

(19) 


(24) 


where the subscript (°°) variables are suitably selected reference values for the dynamic and thermo- 
dynamic variables, L is the characteristic length dimension, and the Lewis and Prandtl numbers are 
to be defined. (Note: in problems corresponding to very slow viscous flow, the pressure, equation 
( 1 8), must be non-dimensionalized by a viscous rather than dynamical reference state. ) 

Substituting equations ( 14) - (24) for the dimensional variables in equations ( 1) - (4), and 
deleting the asterisk notation, yields the governing differential equations in non-dimensional form. 


(pY a ), t = - [pUj Y a --pp- pY^] ;i + S a 
P’t = -lP u il ; i 

(PUi)-t = -lP u i u j + P 5 ij -]b T « 1 ; j + FF p b i 


(25) 

(26) 
(27) 


(pH - (Ec) p), t = - [puj H-^ Tijuj- 


RePr 


MH.il 




Re Pr n 


PQ 


(28) 
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Equations (25) - (28) introduce the important non-dimensional parameters for fluid dynamics as 


Reynolds Number: 

Poo Uqo L 
Re 

Moo 

(29) 

Prandtl Number: 

c 0 M 

(30) 

Eckert Number: 

Ec = U °° 

(31) 


C P T oo 
rOO 


Froude Number: 

USo 

Fr2 u 

(32) 

Lewis Number: 

Iii ir Pr 

(33) 


In equation (27) the assumption that the applicable body force is gravity is evidenced in the term 
containing bj, the unit vector parallel to the gravity vector. In equation (28) Pr,^ is equation (30) 
evaluated at the reference conditions. For a thermodynamically perfect fluid, 


Ec = (7-1)M 2 00 

where 7 is the ratio of specific heats and is the reference Mach Number, defined as 
U M 


(34) 




VtRT c 


(35) 


DIFFERENTIAL EQUATION DEVELOPMENT FOR THE DESIRED 
PROBLEM CLASS IN FLUID MECHANICS 


The solution of equations (25) - (28) represents a formidable task. It is desired to restrict the 
generality of the description to the degree that the concomitant mathematical advantages are appli- 
cable to non-trivial problems. The essence of the development is restriction to two independent 
space dimensions, but in which three-dimensional flows may exist. 

Continuity Equation J 

i 

From mathematical physics it is known that a vector field is completely defined when its 
divergence and curl are known. The vector field of fluid mechanics is velocity, and equation (26) 
defines the divergence of velocity in terms of the density. Identifying the mass flow vector gj as puj, 
the divergence of gj vanishes for either incompressible or steady compressible flows. Henceforth, 
restrict consideration to these two distinct problem classes. 

The vanishing divergence of gj allows specification of gj as the curl of some vector potential 
'Pj. No analytical benefit accrues from this specification (as does occur, for example in Maxwell’s 
equations) in fluid mechanics, due primarily to the nonlinearity of equation (27). But a significant 
numerical benefit can occur if the mass flow vector gj is planar, or if the three scalar components 
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of gj are independent of one independent variable, as occurs in problems possessing axisymmetry. 

In either case, equation (26), for incompressible or steady compressible flow, is identically satisfied 
by the x 3 scalar component of the vector potential 'Fj, in the form 



pui = T «3fl*3;j + PM i3 


(36) 


In equation (36), e 3 jj is the Cartesian alternating tensor, and J is the determinant of the space metric. 


Since all problems of physical interest can be spanned by rectangular, cylindrical and/or 
spherical orthogonal coordinates, equation (36) can be written in terms of a scalar function i//, and 
the gradient operator, as 

pu; — e 3 jj 0, • + pu 3 (37) 

rsin a 0 J 


In equation (37), a is non-zero only for spherical coordinates when it is unity, r is set equal to unity . 
(and u 3 = 0) for rectangular coordinates, and * 3 corresponds to the azimuthal angle (0) in both 
cylindrical and spherical coordinates with domain 0<x 3 <27 r. 


Equation (37) identically satisfies equation (26), in the three Cartesian coordinate systems, 
for both incompressible and steady compressible flow. For example, in spherical coordinates, 


- j 

(puj) . = — — e 3 ij 0,: + pu 3 8 i 3 
;i rsm0 J ;i 


~ (— — 0>j) . TT ^,l) . + (PU 3 ) 

rsirup yrsin <p J ;1 \rsin0 / ;2 


—r- V'.i) + (pu 3 ) 

rsiwp / ;2 »3 J 


Identifying xj as | r, <p, 0 } and transforming the vector derivative to partial differentiation in spherical 
coordinates, yields 


(pUjl-j 

1 


— !_(?_L 

r 2 sin <p \rsin0 


1 r 2 sin0 2 90 ' 

r 90 / > r 


1 rsin0 

rsin0 


£)* + (PU,)># ] 


1 


r 2 sin <p 
= 0 


' 9 2 0 d 2 0 + Q " 

dip 9r dr dip 

l . 


since pu 3 is independent of x 3 , i.e., 0. The proof in cylindrical and rectangular coordinates is similarly 
direct. 
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Compatibility Equation 

With the zero divergence of the mass flow vector gj ensured by identification of the stream- 
function, it remains to obtain its curl. It is useful to identify the vorticity vector 

n i 2 y e ijk u kj ‘ (38) 

The x 3 scalar component (co) of in the three orthogonal coordinate systems, is, 

a) = — e 3 i j u j ;i (^9) 

r® 

with a interpreted as before. 

The vorticity scalar component u> will be employed as an auxiliary dependent variable; there- 
fore, it is necessary that it be related to \p, the primary dependent flow variable. Substituting for uj 
from equation (37), the desired relation is 


co = 


- — ( — 5— *, k 

r a ' prsin a 0 


* 


(40) 


where use has been made of the skew-symmetric properties of the alternating tensor contraction 
e 3 ij ^3jk = ‘ 6 ik 

In rectangular coordinates, equation (40) reduces to the familar Poisson equation involving the 
Laplacian operator. In the other two coordinate systems, equation (40) lacks certain metric co- 
efficients of being the Laplacian operator. 

J 

The Curl of the Navier-Stokes Equation 


Since the problem class is restricted from full-dimensional generality, a specification of the 
x 3 scalar component of the curl of gj yields the desired differential equation. Using equation (27), 
the curl equation is . 


eaki 


P Ui 


i’t 


pUjUj 


P 5 ij 


1 

Re 



_1_ 

Fr 


Pbj 



(42) 


Since the pressure and (gravity) body force fields are conservative, their respective curls vanish iden- 
tically. (The elimination of pressure as a coupled dependent variable is a prime computational fea- 
ture of the transformation to vorticity - streamfunction variables.) Hence, equation (42) can be 
simplified to 




ki * " ‘ 63 ki £ pu i u i ' T u j 3k 


(43) 


ro 



It is required to transform the velocity vector in equation (43) to terms in the x 3 scalar 
components of velocity, vorticity and streamfunction. Considering first the convection term, from 
equation (26), for incompressible and compressible steady flow, obtain 

faki (p u i u j)jk = e 3 ki ^pujUij^ ( 44 ) 

Substitute equation (37) for puj and up Perform the indicated vector differentiation, and using the 
skew-symmetric properties of the alternating tensor contraction 


£3 ki ^3 jC ~ Skj^iC'^kC 5jj 


the convection term becomes 


e 3ki Cp u i u j ) ;jk - e 3 ki 




+ pu 3 e 3 ig 


jl± \ 

prsin a 0 J 


;3J 


(45) 


(46) 


which involves only the dependent flow variables go, ^ and u 3 . Since the vector differentiation has 
been performed where appropriate, the remaining alternating tensors merely restrict the range of sub- 
scripts and provide for correct signs. 

Consider next the term in equation (43) involving the stress tensor Vy. At this juncture, 
assume that either, (1) the fluid is Newtonian, or, (2) that turbulent flow is adequately characterized 
by the laminar flow equations with the molecular viscosity replaced by a scalar turbulent “eddy vis- 
cosity.” In both cases, equation (5) defines the stress tensor constitutive relationship to velocity 
gradient and viscosity. (Note, that since these equation developments begin with the stress tensor 
explicitly, other stress-strain laws could be invoked, as might occur for turbulent flows where the 
eddy viscosity is not adequately represented by a non-directional scalar quantity.) Substituting 
equation (5) into (43) yields 


e 3 ki T ijjk - e 3ki 


P (uij + u j;i) -yM ufi 


# 6 ijj jk 


(47) 


The last term in equation (47) vanishes identically from symmetry considerations. Expanding the 
remaining terms, and noting that the uj-jk term vanishes identically in Euclidean space, the stress 
tensor expression reduces to 


r 


e 3 ki r ij Jk - e 3 ki 


1 


P,k (Ujj + Uj;i) + H Uijk 


(48) 


11 


1 



Considering equation (48) in detail, the last term can be written as . 

P (e s ki u i^c) J j 

Using equation (39), obtain equivalently 


«,ki = (p (r a w),j) -j 


( 49 ) . 


In the remaining terms in equation (48), substitute from equation (37) for the velocities uj and'uj. 
Performing the indicated vector differentiations, and using the symmetry properties for alternating 
tensor contractions, obtain 


£3 ki 


M.k (ui-j + ujj) 


j = -(P^w) 





(50) 


Combining equation (49) with equation (50), and expanding the last term for clarity, the stress 
tensor term in equation (43) becomes, in terms of vorticity and streamfunction, 


‘ ki Tij Jk = ;j 


■+ p,k 63 ki 5 js u 3 ;i 


i 

- ;j 


(51) 


There remains the left-hand side of equation (43) which is non-vanishing, within the assump- 
tions regarding the validity of 1 p, only for incompressible flow. For this case, 


P e 3ki u i;kt - (52) 

Combining equations (46), (51) and (52), the x 3 scalar component of the curl of the Navier- 
Stokes equation (27), written entirely in terms of the vorticity, streamfunction and x 3 component 
of velocity u 3 , becomes 


(pr a co), t 



2 


2 



1 r 

+ — 1 
Re 


.+P, k e 3ki 5j 3 u 3;i (53) 

J \ prsin a $ / ;j J ; j 


The form of equation (53) is desired for computational purposes, since the velocity vector and 
pressure have been explicitly eliminated. However, conceptual clarity results from transformation of 
the terms containing t//, k into velocities. Contracting equation (37) with e 3 j k and using equation (41) 
yields - 

** 


prsm a <t> 


= fsjkuj 


(54J 
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Removing the i^, k the terms in equation (53), the equivalent form with velocities is 
pr a co, t = - |\pu k r a co) ;k +~ e 3 ki( u £ u £);i P>k + e 3ki(pMj 3 ( u r u 3 5 i3) ; j 

+ ' P>k ' 2 P’j e 3£j ;k J ;k 

The convention on x 3 in each coordinate system allows simplification of terms involving pu 3 . 
For rectangular, cylindrical, and spherical coordinates, respectively, define 

Xj= {x,y,0[ 

= l z - r > 0 l ( 56 ) 



Then, for example, in cylindrical and spherical coordinates, respectively obtain 


e 3ki (P u 3 6j 3 ( u j ' u 3 ^ i3 ) ;j );k 


fjfc 


(pu 3 ) 2 

>z 



(57) 


The x 3 Component of the Navier-Stokes Equation 

The determination of u 3 , in terms of co and \p,, is obtained from solution of the x 3 component 
of the Navier-Stokes equation, which is 

(pu 3 ), t = - (pu^a)^ - p,3 + T 3 j\j (58). 


assuming the body force distribution has no x 3 component. The x 3 derivative of pressure is included 
in equation (58) to account for flow processes through “corkscrew” devices.- 

Using equation (37) to replace pu^ and proceeding through the details of the vector derivative 
of u 3 , equation (58) becomes 


(P u 3 )> t = -g 3 ij ( ~ ~~4 U 3 ) -/>U 3 U 3;3 - p, 3 

\ rsin a 0 / 

+ Ri [(P U 3^);k+y.PU j;j 3 + M, k U k;3 ] 

Note that uj ;3 does not vanish identically. For example, in cylindrical coordinates 


(59) 


1 

U 3 ; 3 = — Uj (60) 


1 

u 3 ; 3 3 ~ -~r U 3 
r 


(61) 
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For incompressible flow, the transient term may be retained; otherwise, only the steady state solu- 
tion is correctly specified. 

The Species Continuity Equation 

For multi-component flow, equation (25) describes each species mass fraction with respect 
to mutual convective, diffusive and reactive processes. Transformation to the desired computational 
form simply involves replacement of puj by the streamfunction definition, equation (37). Since the 
divergence of equation (37) vanishes, and Y a is not a function of x 3 , the desired equation becomes 



The number of equations (62) equals the total number of species present, and the transient solution 
is correctly specified only for species of identical molecular weight. 


The Energy Equation 


The form selected to express energy conservation, equation (28), is particularly well suited 
for numerical computation in the problem class being considered. The explicit pressure influence is 
contained solely within the temporal dependence. Since the pressure, as an explicitly coupled de- 
pendent variable, has been eliminated from the momentum equations as well, solutions to steady- 
state problems for either incompressible or compressible flow are possible, without solution fox the 
pressure field, provided that the implicit influence of pressure as a thermodynamic variable can be 
suitably approximated. This advantage is particularly noteworthy, since most numerical algorithms 
for fluid mechanics, which contain pressure as an explicitly coupled dependent variable, are specially 
designed to “half-step” iterate between the flow field and the pressure distribution. Inefficiency or 
numerical instability characterizes the attempt to march both variable types simultaneously. 


To cast equation (28) into the desired computational form, replace the velocity vector by the 
streamfunction definition, equation (37). The convection term becomes 



g 3ij 

rsin a 0 


+ pu 3 



(63) 


Using equations (5) and (6), the irreversible work term is 


Tjj Uj = n 


1 2 
y (UjUj);i + UjU i;j ■ y U i U k;k 


( 64 ) 


Using equations (37) and (41), the first right hand term becomes 



Using equation (45) also, the second term becomes 


u j u ij 


1 

T 


*’k 




*4 


v prsin a $/ ;i prsin a 0 \prsin a 0/ 3 
jk ^ »k 


| N 

■ — (Uj )* - 3 + U 3 


w . n , 

prsin<0 . 


( 66 ) 


The final term in equation (64) can be transformed using equation (26), to 

u i u k;k = - u i u k 8n ^*k < 67 > 

Its desired computational equivalent is 

uu - 1 / tl i \ V 1 t 1 / ^’k ^>i \ / 1 \ 2 e 3kg^»g /I \ 

1 k,k 2 \rsin a 0/\p /.j 2 \(rsin a 0) 2 / \p / , k rsin a <£ \pj 5 i 3 (68) 


Combining equations (64) - (67), the irreversible work contribution to the energy equation 
becomes 
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The computational form for the energy equation thus becomes 


0>H), t = (Ec) p, t 


/ C 3 kj *.j 
\ rsin a 0 


+ PU 3 



H 


Re 



pEc 

/ *>i ' 

’ 

*’k 

( y 


"rT 

yprsin**^ 

/ *' 

prsin a 0 

\prsin a 0/ ;j 


pEc 

iRe 

u 3 ;k + 

( u » ;» 

+ 2u 3 j 

8 k ,' 

prsin a 0 J 



For compressible flow, the time derivative terms must be set to zero. Equation (70) is valid as 
written for transient incompressible flow. 


Recovery of Pressure 


The pressure distribution is recovered from the vorticity-streamfunction characterization of 
this problem class by enforcing linear momentum conservation. Equation (27) can be transformed ' 
to the Laplacian on pressure by an additional differentiation, and solved as boundary value problem. 
Since this is not a well-posed description for pressure, the preferable approach may be to integrate 
equation (27) after contraction with the infinitesimal displacement vector, dxj, yielding 



-n 


PU ‘ U > -Re 


4 






(71) 


Note that the left side of equation (71) is the integral of a perfect differential and thus independent 
of path. Hence, the pressure at any point in the field can be determined, in comparison to some re- 
ference value, by integration over an arbitrary (the easiest) path between the two points. Thus, the 
entire pressure field need not be computed to obtain selected pressures, for example, at an exterior 
wall. Upon noting that the integrand in the last term of equation (71) is Xj independent, and denoting 
integrals of perfect differentials as A, equation (71) becomes 



;j dx i + k / T y ;J dx i ' Ji / pUidXi + "Fr pAXb 


(72) 
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To consider the convection contribution to equation (72), replace pujuj by equation (37). 
Performing the indicated differentiation, using equations (41) and (45), and noting x 3 independence 
of some integrands, obtain 


/ (^ u i u j) j dx i = A 


i / ** V 

-[ 

1 

— i 

P Irsin a <pJ 

J 

p(rsin° £ 0) 2 j 


dX: 


/ u 3; 


• e 3 i£ 'P >£ 


3 dX: 

rsin a <t> 


+ Ax 3 


„U , , t L 


e 3j \ 
rsin a <£ / ;j 


(73) 


Recall that u 3 ;3 is transformed to equivalent terms in i p for each particular coordinate system. 

The same procedures are used to transform the stress tensor contribution to computational 
variables. Integrating selected terms by parts, the following form can be obtained. 



dX: 


A 



g -»jk.*.k 

prsin a 0 


;j 



e 3jk^»k 

prsin a 0 


Snp ,j 


y 



e 3jkV/, k 

■r*»i — — 

prsin®^ 




e 3 ik ^ ’k 
prsin a tf> 


;j 


dx i 


Ax 3 


(" u » j) 


/ 


Mtl 3 ;j'i dXj 


(74) 


Substituting equation (37) for the integrand of the time dependent contribution to equation (72), 
and combining terms, the following form results for determination of the pressure at any point in 
the field. 
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Ap = A 


1 

7^3jk^-k\ 

Re 

_\ prsin “0 / 


2 e 


3 jk^»k /J_\ 
in a <t> \ P / , j _ 


3 rsin a 0 \ p 


M—, j +' jT 

p \rsin a 0 / ^r 


M / ^ e 3 ik ^ »k \ 

Re \ prsin a 0/3 


- ^>i 


e 3jk^»k 

prsin a 0 


px 5 


P,j 


C 3jk^.k 

prsin a 0 


*>i tA 
2 

p(rsin a ^ 


;j dx i 


►/ 


1 es 

P u 3 ;3 i - u 3;3 


3 ik 0 »k _ _3_ / e 3 jk ^ >k\ 

rsin a 0 3t y rsin a <t> ) 


dx i 


+ Ax 3 


^ (> u »a) ;j - ^ u % ;3 1 ( u 3 . 

Re 7 \ rsin “0 / ;j 


(75) 


In equation (75), the first curly bracket contains terms which are point dependent, the result of 
integrating perfect differentials. The second bracket requires integration, and hence selection of 
path. The third bracket contains terms which require evaluation dependent upon the selected 
coordinate system. The last bracket merely requires evaluation, since all variables are x 3 -independent. 


FINITE ELEMENT SOLUTION ALGORITHM FOR THE 
NAVIER-STOKES EQUATIONS 

The desired form for the differential equation system governing the mechanics and thermo- 
dynamics of a viscous, variable density fluid in two independent space coordinates has been estab- 
lished. The system of elliptic partial differential equations, which may be initial value as well under 
certain restrictions, requiring solution includes equations (40), (43), (59), (62), and (70) for stream- 
function ( 1 p), vorticity (to), swirl velocity component (u 3 ), mass fraction (Y a ), and stagnation en- 
thalpy (H). Equation (75) is an algebraic equation for determination of pressure (p) in terms of the 
flow computational variables, and the density (p) is established via equation (7). Closure of this 
complex equation system is obtained by specification of the thermophysical and transportive pro- 
perties of the fluid, that is, laminar flow molecular viscosity, thermal conductivity, and specific heat. 
The formulation is equally valid for steady turbulent flows, described by the time-averaged Navier- 
Stokes equations, by establishing an “eddy viscosity” model based upon the kinematics of the flow 
field. The solution algorithm to be established assumes a generalized viscosity description that may be 
an explicitly prescribed function of dependent as well as independent variables. 
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The uniformity of the derived differential equation description for the computational de- 
pendent variables is evident as specific variants of the nonlinear elliptic partial differential equation 
of mathematical physics. Identifying q as a generalized dependent variable, the uniform description 
belonging to each of the equations (40), (43), (59), (62) and (70) is 

N(q) = K[Kq, k ] ; k - [e 3ki ■ q] ; k + ftq.q.j^Xj) - g(p,q,t) = 0 (76) 

rsin“0 

In equation (76), N(q) is the specified differential operator for q, the first right side term is the 
elliptic operator, the second term is the convection operator, and f and g are specified functions of 
their arguments, see Table 1 , corresponding to q identified with each dependent variable. Since 
equation (76) is elliptic, unique solutions are obtained only after specification of boundary con- 
ditions, the most useful and general form of which relates the function and its normal derivative 
everywhere on the closure, 3R, of the solution domain, R, as 

rj(q) = aO ) (Xj,t) q (Xj,t) + a( 2 ) (xj.t) q (xj,t) , k n k - a( 3 )(Xj,t) = 0 (77) 

In equation (77), the a(*) (xj,t) are user-specified coefficients enforcing the boundary conditions 
which may be time dependent, the superscript bar notation constrains Xj to 3R, and nfc is the local 
outward pointing unit normal vector. For g nonvanishing in equation (76), or for use as an esti- 
mated distribution to start an iteration to steady state, an initial distribution is also required; hence, 
assume given throughout R 

q(Xj,0) = Q 0 (x i ) (78) 

Formation of the finite element solution for the two-dimensional Navier-Stokes equations 
is thus obtained by establishing the algorithm for the equation system (76) - (78) identified for each 
dependent variable requiring solution. The theoretical foundation is the method of weighted resi- 
duals (MWR) (Reference 1) applied on a local basis. Discussion on the formulational aspects of this 
algorithm are reported for a considerably restricted problem class within the Navier-Stokes equations 
(References 6, 8) as well as for three-dimensional boundary-layer-type flow fields (References 7, 10). 
Since equation (76) is valid throughout R, it is valid in disjoint interior subdomains Rm> described 
by (xj,t)e R m x [0, °°) called “finite elements,” wherein UR m = R. Form an approximate solution 
to equation (76) within R m x [0, »), called q^ (xj,t), by expansion into a series solution of the form 

q m<V> s |, t L* <*i> 0? <‘> = iUXi)j T jQ(t)) m (79) 

wherein the functionals L k (xj) are members of a function set that is complete in R m , and the un- 
known expansion coefficients, Q k (t), represent the time-dependent values of q^j (Xj,t) at specific 
points called “nodes” interior to R m and on the closure, 3R m . Equation (79) is a scalar, the second 
form introduces the column matrix notation and its transpose (superscript T), and the method of 
selection of the L k is user specifiable and problem class dependent. To establish the values taken 
by the expansion coefficients in equation (79), require that the local error in the approximate solu : 
tion to both the differential equation (77), N(qm), and the boundary condition statement, T?(qm)> 
if applicable, that is, 3R m coincident with 3R, equation (78), be orthogonal to the functional set in 
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, TABLE 1 

COEFFICIENTS IN GENERALIZED, ELLIPTIC. DIFFERENTIAL 
; EQUATION DESCRIPTION, EQ. (76) 



rsm“0 






equation (79). This is accomplished by enforcing the local variant of the Galerkin criterion of 
classical MWR, that is, premultiplication of equations (77) and (78) by the approximation functionals, 
L k (xj), and integration over the respective domains. Employing a Lagrange multiplier, X, these 
equation sets can be combined to yield 

/ {Lix^m^XlT - \j> (LfxpItfq^Wo a 0 . (80) 

R m dK m 

The number of equations (80) is identical to the number of node points associated with the finite 
element, R m , that is, the number of elements in the column matrix { Q(t) ) m . 


Equation (80) forms the basic operation of the finite element solution procedure. Establish- 
ment of the global solution algorithm, and determination of X, is accomplished by evaluating equa- 
tion (80) in each of the m finite elements of the discretized solution domain, R, and assembly of 
these Mxn equations into a global matrix system by Boolean algebra. The rank of the global system 
is less than Mxn by connectivity of the finite element domains as well as by enforcement of those 
boundary condition statements on 9R where a( 2 ), equation (77), vanishes identically. Determination 
of X is accomplished by observing that the lead term in equation (76), when incorporated into equa- 
tion (80), can be integrated by parts; alternatively, using the generalized divergence theorem (Appen- 
dix A), obtain 

/ jL( Xi )| K [Kq^ ] dr = k $ | L(xj) ) Kq^ n k da' - k/ j L(xj) ) , k Kq^ dr 
*m ** dK m K m 

' (81) 

On those portions of the closure 9R coincident with 9R m , equation (81), the corresponding segment 
of the closed surface integral will cancel the boundary condition contribution, equation (80), by 
identifying X with k of equation (76) and a( 2 ) with the generalized diffusion coefficient, K. From 
Table 1 , for the Navier-Stokes equations, k typically involves the Reynolds number and a( 2 ) is the 
viscosity, p, except for the vor.ticity-streamfunction compatibility equation (40), where they are 
unity and inverse density. The contributions to the closed surface integral, equation (81), where 
dR m is not coincident with 9R remain for evaluation. Hence, combining equations (77) - (81), the 
globally assembled finite element solution algorithm for the. representative elliptic partial differential 
equation system description for the Navier-Stokes equations becomes 



Using Boolean algebra as the assembly (summation) operator, the rank of the global matrix 
equation system (82) is identical to the total number of node points within the solution domain R, 
and on the closure 9R, at which the dependent variable requires solution. Equation (82) is either 
a first-order ordinary differential, or algebraic equation system dependent upon whether g, equation 
(76), vanishes identically. In either event, the equation system is large order and the matrix structure 


I 
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is sparse and banded about the main diagonal, with bandwidth a function of both selected discreti- 
zation and the order of the employed approximation function, equation (79). Solution of the 
algebraic system can be obtained by many procedures including direct inverse, Gauss elimination, and 
equation solver techniques, for example, banded Cholesky. If the equation description is ordinary 
differential, solution is obtained using either an explicit or implicit numerical integration algorithm 
for an initial value problem. 


FINITE ELEMENT MATRIX GENERATION 

Implementation of the finite element theory into a computer code involves assembly of 
Eq. (82) for the required dependent variables, hence evaluation within each element of the discre- 
tized solution domain. The vast experience in structural mechanics has proven the viability of trun- 
cated power series as approximation functions. Highly tailored finite elements have been established 
for specific problem classes in elasticity. But for now, the use of finer discretizations using linear 
(simplex) approximation functionals, Eq. (79), appears'useful for the broader problem class. The 
intrinsic finite element shapes for one-, two, and three-dimensional spaces, spanned by simplex 
functionals, are the line, triangle, and tetrahedron. The preferred orientations for these elements 
in a local coordinate system are shown in Figure 1 . 

Accurate determination of the element matrices of Eq. (82) is mandatory, and involves evalu- 
ation of various-order moment distributions over the domain of the finite element. A straightforward 
approach involving numerical matrix multiplications, after definition of a transformation matrix 
(Reference 8), is generally adequate for linear equations; however, it has proven highly susceptible to 
accumulated round-off for the explicitly nonlinear matrices common to the Navier-Stokes equations. 
An alternate approach has been developed which expands upon the natural coordinate function 
description for finite element solution of the energy conservation equation (Reference 1 1). These 
functionals, which are an adaptation of the area coordinates of structural mechanics (Reference 
1 2) are a linearly dependent set of normalized functions that are orthogonal to the respective closure 
segments of the finite element domain. For an n-dimensional space, there are n+1 simplex natural 
coordinate functions. Table 2 contains the implicit definition of these functions in their respective 
spaces. (The three-dimensional descriptions, while inappropriate for the particular Navier-Stokes 
equation system being studied, are included to complete the finite element theoretical foundation 
for future reference.) As the direct consequence of the definition, the natural coordinate functions 
vanish at all node points of the finite element except one where the value is unity. Hence, these 
functions are the elements of the approximation function column matrix, j L } , Eq. (79). Of 
particular impact for solution of the Navier-Stokes equations is that integration of arbitrary -order 
products of scalar components of the j L } , over the domain of the finite element, are analytically 
evaluable in terms of the exponent distribution, Table 3, the determinant of the equation system 
defining the Lj, Table 2, and the dimensionality of the space of the problem. 

It remains, to evaluate the solution algorithm, Eq. (82), for the Navier-Stokes equations using 
the natural coordinate function description, which involves matrix formation over one- and two- 
dimensional finite elements. (The one-dimensional element is needed to evaluate the domain closure 
integrals of the solution algorithm.) The evaluations are presented for the two-dimensional flow of a 
fluid with arbitrary density and viscosity distribution, and the axisymmetric flow of a constant prop- 
erty fluid. Additional integrations by parts have proven useful for the functions f(q), Eq. (76), and 
the mathematical consequences of these operations are established in Appendix B. 
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Figure la. One Dimensional Space - Figure lb. Two-Dimensional Space 


*3 




Figure lc. Three-Dimensional Space 


Figure 1 . Intrinsic Finite Element Domains for Simplex Approximation Functions 
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TABLE 2 

IMPLICIT DEFINITION OF SIMPLEX NATURAL 
COORDINATE FUNCTIONS 



TABLE 3 

INTEGRALS OF NATURAL COORDINATE FUNCTION 
PRODUCTS OVER FINITE ELEMENT DOMAINS 


Dimensions 



Integrals* 


D 


n, ! n 2 ! 
(n+ n ] + n 2 )i 


2 


3 



n 2 n 3 
L 2 L 3 dxdy 


n ? n 3 n 4 
L 2 L 3 L 4 dxdydz 


= D 


n i ■ n 2 I n 3 ! 
(n+nj +n 2 + n 3 )! 


rij I n 2 ! n 3 ! n 4 ! 

(n+ nj +n 2 +n 3 +n 4 )! 


D = Determinant of coefficient matrix defining the 
natural coordinate system, see Table 2, 

n = Dimensionality of the finite element space 
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Planar Finite Elements for the Two-Dimensional Navier-Stokes Equations 

With the advent of the natural coordinate function description for the Finite element approxi- 
mation functions, matrix evaluation is straightforward. Moment generation is invariant under both 
coordinate translation and rotation in Euclidean space spanned by a rectangular Cartesian basis, 
(Appendix C). All computations are formed in the local (primed) coordinate system, Figure 1 , 
defined by the tensor transformation law 

x 'j = a jj x j + r j (83) 

where r, is the position vector to the origin of the primed coordinate system, and the ay are the 
direction cosines of the coordinate transformation. The integration kernels for two-dimensional 
space, Eq. (82), are 

dr = t m dxi dxi . (84a) da = t m dx', (84b) 

where t m is the thickness of the m^ 1 finite element. (Variable element thickness allows solution of 
predominantly two-dimensional flows in irregular depth channels, for example, rivers.) Since no 
difference exists between differentiation of scalars and vectors, the semi-colon operator reduces to 
the comma. Furthermore, no u 3 velocity component can.exist in rectangular Cartesian space. 

All expressions but the third, Eq. (82), are standard for all dependent variables. For the 
generalized diffusion term, using Eq. (79) for both and the generalized diffusion term, K m , and 
transposing scalar quantities, obtain 

* f R ld ■k K "'‘lm, k dr " “I' M mW hl, k jL)T k lQl m dr 

= «|K|£ |B!0| [B21IS]|Q| m (85 ) 

The standard matrices, as used in Eq. (85) and the following, are defined in Table 4. If j Q } in Eq. 
(85) is identified, for example, with vorticity, then k - 1 equals the Reynolds number, Re, and the 

elements of | K } m are the element node point values of viscosity, M- The second term of Eq. (82), 
which accounts for convection processes, becomes 

f Id e 3 k j 7~a7 ; k d T = f |'/'| m I ’i e 3 ki IM > k |Q| 

= |BI°I !*!■„ [B211A] |Q| m (86) 

Equation (86) is eligible for integration by parts; however, no arbitrariness exists for its evaluation, 
Appendix B, so no benefit accrues from that operation. In contra-distinction to this, the contribu- 
tions to Eq. (82) stemming from the final term, involving a closed integration on 0R m , are non- 
vanishing for all R m , but can be made to cancel in pairs upon Boolean assembly of the global algor- 
ithm, see Appendix B. Hence, the sole surface integrals requiring evaluation occur when 3R m and 3R 
are coincident, and a^ 1 ' and/or a m ^ 3 are non-vanishing. Therefore, the only contribution to Eq. 
(82), from all surface integrals, reduces to two terms, 
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TABLE 4 

STANDARD FINITE ELEMENT MATRIX FORMS FOR SIMPLEX FUNCTIONALS IN 

TWO-DIMENSIONAL SPACE 


MATRIX 

nameO 


(B200S1 


[B21IS] 


MATRIX 

FUNCTION 


(L| {L}T dT 


IM. k {U >k 


MATRIX EVALUATION* 2 * 


2 1 ll O) 


far [' ■ :J 

/X2P3 V 2 X2P3 /X2P3 \ / X2P3 \ 

\XIP2 */ ’ XIP2 \X1P2 "/’\XIP2 / 

J_1 \ 2 /X2P3 \ 2 / X2P3 \ 

\X2P3 / ’ ( XIP2 / \ X I P2/ 


IB2IIA] e 3ki {L|,jL}7 k 


{L{ dr 


[A200S] 


W |L} T do 


% 2 l 0 
2 0 
0 


da 


(1) Matrix names are a 6 digit code covering dimensionality, nonlinearity, degree of 
differentiation and special matrix properties, as [a, b, c, d, e, fj where: 

a = A, B, C for spaces of one-, two-, and three-dimensions, 

b = number of natural coordinate functions appearing in integral or matrix, 

c, d, e = (0, 1) Boolean counters indicating (no, yes) differentiation of each function, 

e or f = S, A, A for matrix symmetric, antisymmetric or general. 

(2) A m = l A (XI P2) (X2P3), the plane area of the (triangular) finite element, where the 

four digit codes signify node point coordinates in the primed system, for example, 

X1P2 = the x, prime coordinate of node 2, 

X2P3 « the x 2 prime coordinate of node 3, 

t m = average thickness of the finite element, 

C m = length of sink along which the boundary condition is 
applied (=X1P2) 


(3) Symmetric matrices are written in upper triangular form. 
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(87) 


«/, R |L| a m O q > = «a m O)[A200S]|Q| m 

l L l a m <3) da = Ka m ( J) | A 1°1 (88) 

where the a m (^ are user-supplied boundary condition constraints. 

One more standard matrix exists should the analysis be performed for the transient motion of 
an incompressible fluid, that is, hydrodynamics. The function g^, Eq. (82), is non-vanishing for all 
dependent variables excepting streamfunction. The numerical evaluation for vorticity and mass 
fraction takes the form: 

L Id 8^dr = / |t) p m |L| T $|CXt)| m dr = p m [B200S] |Q(t)| ' m (89) 

m ^m 

In Eq. (89), the superscript bar signifies an element-averaged value for density, while the prime 
denotes that the dependent variable matrix contains the time derivatives of the node point 
values of the respective functions. Hence, in this case, Eq. (82) is a system of ordinary differ- 
ential equations. For the energy equation, a more convenient form exists in terms of temper- 
ature as the time dependent variable. Hence, obtain 



= Pm c P IB200S| hOlm + | £ -| B l 0 l 1*1 „ IB21.1S] (90) 

Since the mechanics and thermodynamics of incompressible flow are separable, {</(} ' m can be in- 
dependently evaluated using a finite difference formula for the derivative of a tabulated function. 

Thus, the second term in Eq. (90) is at most a source term, usually of negligible magnitude, since 
the Eckert number is typically small for the more common incompressible flow configurations. 

Several additional matrix evaluations are required for fm to complete the description of 
Eq. (82) for specific dependent variables. Taking them in the order of occurrence in Table 1 (f^j is 
the approximation to f in R m ), for streamfunction, the second term simply cancels the convection 
operator of the standard equation form, Eq. (82). For the remaining term, obtain 

L {Ljo^dT- = [B200S] \n\ m ( 91 ) 

For solution of vorticity, the terms involving u 3 are discarded, and the remaining terms can be 
effectively integrated by parts to render the forms more tractable. Considering the density term, con- 
tributions stemming from the resultant surface integral can be made to vanish in pairs bn all dR m inte- 
rior to 3R upon assembly of the global algorithm, as discussed in Appendix B. For locations where 
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3R m coincides with 3R, no arbitrariness exists for evaluation of the surface integral. However, 3R 
coincides with \p equal to a constant everywhere except at mass flow injection points. At these 
locations, the density may usually be assumed constant, and either condition is sufficient to render 
the matrix identically zero. So, the sole term requiring evaluation is 

L l*L[B2HA]|-J r ) In (92) 

K m X P z 

where p -1 is selected as the preferred primative computational variable due to its universal appearance. 
Thus, the elements of j-Ll are the inverse of the node point densities of R m . 

After integration by parts of the viscosity gradient terms, no surface integral evaluations on 
interior dRm are required, based on the analysis used for density. For those segments of 3R m , 
coincident with 3R and for no mass injection, the flow is laminar, hence viscosity is only weakly 
temperature dependent. At injection points, the density may be assumed constant, and either con- 
dition renders the surface integral terms identically zero. So, for variable viscosity, the sole contri- 
buting terms from fm for vorticity become 

T^/ Rm W. k ^k< + 2^jC, j (7).k)^ IB21 ]S] |«| m 

*TT [B211S1 {£}; m |m| ^[B211SJ |*) m (93) 

where the elements of the column matrix |jli} are the node point values of viscosity. 

The sole term in f^, for the species continuity equation is that due to sources of mass. This 
term could be conveniently expressed as either element or node point dependent. The respective 
matrix forms are, 

/ |L} S a dr = S® <B10} 

R m 

= [B200S] |S a | m 

For energy conservation, Eq. (70), the entire viscous dissipation term may be integrated by 
parts, and the resultant surface integrals on interior 3R m vanish in pairs upon assembly and rendered 
zero on exterior surfaces. The various terms are combined to yield the following form 



= ( t H ra (B211S]{i} m |^|^B200S)(I} n , 

- l*l^[B211Sl{i} m [B211S] |*| m |M|^[B200S](I} m ) (95) 

/ 


(94a) 

(94b) 
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Finally, the source term, pQ, in the energy equation would typically be node point dependent; 
hence 


/ |L}pQ m dr = [B200S] {pQ} m 


(96) 


It remains to establish the boundary condition specifications, a^, Eq. (77), for the 
derived dependent variables, streamfunction and vorticity, in terms of known velocity and 
density distributions. The practical engineering specification is velocity scalar components per- 
pendicular and parallel to the domain closure 3R, and mass flux across appropriate closure 
segments. The closure specifically need not be, and oftentimes is not, parallel to coordinate 
surfaces of the global reference frame. Establishment is straightforward; contracting Eq. (37) 
with an alternator, and using Eq. (41), obtain 

~ ^’k ~ e 3ik u i (97) 

Referring to Figure lb, identify the surface between nodes 1 and 2 of R m as coinciding with 
a segment, 3R m , of 3R whereupon a specified velocity and mass flux distribution exists. Assume 
the velocity scalar components perpendicular (uj_) and parallel (uy) to the line are assigned at 
each node point. Contract Eq. (97) with an infinitesimal displacement vector tangent to 3R m 
and integrate to obtain 

(2) (2) (2) 

x i j i x i / x i t 

J ~ ^,, dx i = / „ e 3il u i dx i = j u i dx i (98) 

0 0 o' 


Using Eq. (79) for all variables, the integrals are directly evaluable. Table 3; hence for a spec- 
ified normal velocity on 3R m , the equivalent streamfunction constraint is 



(99) 


Equation (99), in combination with the normal mass flux specification, determines both \p and 

— on 3R m . In Eq. (99), u^ is defined positive for mass efflux from R m , and £ m is the length 
of the closure segment between nodes 1 and 2. 


The counterpart evaluation of Eq. (97), contracted with the outward pointing unit normal 
vector, n^, becomes 


1 1 I 

-*’k n k = C 3ik U i n k * - U || 


( 100 ) 
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where uy is defined positive for alignment parallel to the x',axis, Figure lb 
the equivalent stream function boundary condition constraint becomes 

( 3 ) 

a m “ ' u ll 


The vorticity equivalent of the no-slip, wall boundary condition requires establishment. 
Referring to Figure 2, finite element discretization of R using linear functionals involves 
approximation of arbitrarily curved surfaces by piecewise continuous chords. These surface 
approximations, by definition, are lines of constant streamfunction, even at mass injection 
points, see Figure 2. At each node point on 3R, an approximation to the local normal to 
the surface can be constructed by bisecting the included angle, as shown. Along each of 
these lines the streamfunction-vorticity compatability equation, Eq. (40), locally approximates an 
ordinary differential equation in the direction normal to the surface. After integration by parts, 
the equation takes the form 


. Using Eq. (77), 


( 101 ) 


£ T) 

|ds Jtu(r?)di 7 = 
o o 



( 102 ) 


Using Eq. (79) to represent the functional dependence for all dependent variables, the values 
of vorticity, J2 W , at node points lying on the solution domain closure where no-slip is enforced, 
is expressible in terms of wall and next-interior node point values of the variables as 



(103) 


where £ is the normal displacement of the interior node from 3R. .The vorticity vanishes 
identically on portions of the solution domain closure that are planes of symmetry. For 
closure segments where unknown (outflow) velocity profiles exist, use of the vanishing normal 
gradient of vorticity is appropriate. 

l 

■ . ' ' I 

Equations (83) - (103) complete the finite element solution algorithm specification for 
numerical solution of the steady-state, two-dimensional, Navier-Stokes equations for a fluid 
with variable density and viscosity. For flows where the density is time invariant, the transient solu- 
tion is additionally specified. Should violation of the stated assumptions regarding density and/or 
viscosity variations on 3R occur for a specific problem, the discarded surface integrals can be readily 
evaluated and their contribution to the solution inserted appropriately. 


30 



Figure 2. Establishment of Vorticity Boundary Condition Statement 

Planar Ring Finite Elements for the Axisymmetric 
Incompressible Navier-Stokes Equations 

In axisymmetric space, moment generation is invariant only under coordinate rotation, see 
Appendix C. The preferred element orientation remains in the primed coordinate systems, Figure 1, 
as defined by Eq. (83). The integration kernels become, "respectively, 

dr = 27rx 2 dx 1 ' dx 2 (104a) 

do = 27rx 2 dxj. (104b) 

where x 2 is the radial coordinate in the global reference frame of a point in R m . Before integrations 
can be performed, it is necessary to express x 2 in the primed system; using Table 2 obtain 

x 2 = xP ) + x-Lj (105) 

where x/ * ) is the radial coordinate of node point 1 in the global reference frame. 

The finite element counterpart of the ‘generalized diffusion term, Eq. (82), now requires addi- 
tional consideration for streamfunction. Since the diffusion coefficient contains a metric coefficient, 
Table 1, the resultant differential operator is not the Laplacian in cylindrical coordinates. The explicit 
appearance of a Laplacian can be extracted by definition of a new variable, 0 = 0/X2, whereupon the 
compatibility equation, Eq. (40) takes the form 

-o)= [^0, k ]; k + 0 (p3Tj), 2 (106) 

For the function 0, as well as the remainder of the dependent variables, the generalized diffusion term 
for constant density and viscosity takes the form 
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( 107 ) 


«/ |L|. k K m q^;kf - «/ |L| . k K m |L|T k |Q| m dr 

= »<K m { R } ^ |B10| IB211S] (Q } m 

where K m is an element scalar stemming from the assumed constant properties, the elements of { R } m 
are the node radial coordinates in the global reference frame, and in axisymmetric space, the element 
thickness t m defined in j B 1 0 } , Table 4, is replaced by 27r. 


The convection matrix becomes 



= I BIO } * [B211A] |Q } m . (108) 

which is identical to the two-dimensional form with t m replaced by 2ir. The boundary condition 
matrices for the dependent variables are 


kJ 

J 3 R m' 

1 L ) a m (1) 9m d0 

= Ka m O) 

[A200SA] 

io} m 

(109) 

K / 1 
J dR 1 
m 

L } a m (3) da 

= Ka m ( 3 > 

{a.ioa} 


(110) 


In Eqs. (109) - (1 T0)J the suffix A in the matrix name refers to the axisymmetric counterpart of the 
standard form, see Table 5. In each case, the first term is identical to the two-dimensional form with 
t m replaced by 2tt. 


The transient algorithm is correctly specified for constant density flows. The axisymmetric 
finite element evaluation for g^, Eq. (82), becomes 



( l L Ml-l T -5 r lQl m <lr 



= p m [B200SA] (Qtt)}^ 


(111) 


The added complexity of the [B200SA] standard matrix, Table 5, stems from lack of translational 
invariance. 

' The sole source terms requiring evaluation, for constant property axisymmetric flows are in 
the modified streamfunction equation, Eq. (106). For the vorticlty term, obtain 


f { L | uj^dr = [B200SA] j«} m 

J r> 


( 112 ) 
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TABLE 5 

ADDITIONAL STANDARD MATRIX FORMS FOR SIMPLEX FUNCTIONALS 

IN AXISYMMETRIC SPACE 


MATRIX 

NAME 


MATRIX 

FUNCTION 


MATRIX EVALUATION^ 1 )» 


[A200SA] 


|aioa| 


[B200SA] 


TB200SCJ 


/ iLi(L!T d0 

3R m 


/ fL|do 
9R m 


/ W(L| T dr 

R m 


R 


m 


2irZ m 

“6“l r l 


2 1 O' 
2 0 
0 


r 2’ r l 


1 1 0 
3 0 
0 


27rC m 


r, 1 + 


27rA 


m 


12 r I 


/ {l}{l[Lt^ 3 ) 


2 1 f 
2 1 
2 


^ 111 ) 

r 2~ r i r 

5 L 


2 2 1' 
6 2 
2 


r 3' r l '" 2 1 2 


2 2 
6 


27rA 


m 


12 


’ 2r i + r 2 + r 3 

r l + 2r 2 + r 3 

r l + r 2 + 2r 3 ] 


( 1 ) The lower case r refers to node radial coordinates in the global reference frame. 

(2) Symmetric matrices are written in upper triangular form. 

(3) fB200SCj is a diagonal matrix. 


For constant density, the form for the remaining term is 


/ WW T i»U^d,. 

D 2 

K m 

■ 1820031 ^ I- 

where [B200S] is the rectangular Cartesian matrix defined in Table 4 with t m replaced by 2tt. The 
boundary condition statements for both streamfunction and vorticity are identical to the two dimen- 
sional specifications, with the addition of respective inverse node radii to the denominator of 
Eq. (99) and the numerator of Eq. (103). 



COMOC COMPUTER PROGRAM 

The COMOC (Computational Continuum Mechanics) general-purpose computer program 
system embodies the derived finite element solution algorithm for the general initial value, elliptic 
boundary value equation system, and is coded entirely in terms of generalized non-dimensional 
independent and dependent variables. In addition to producing the results to be discussed, it has 
been'exercised for problems in transient heat conduction (Reference 1 1 ) and the three-dimensional 
boundary layer equations (References 7 and 9). The program consists of four basic Modules. 

In the INPUT Module, the desired discretization is formed by specification of the plane coordinates 
of the triad of vertex nodes for each finite element. The non-vanishing boundary condition 
coefficients a^ , Equation (82), are included as element information for each dependent variable 
to be solved, since in general the different Q’s have individual specifications. The next two Modules 
are basically DO loops on the finite elements of the discretization. In the GEOMETRY Module; 
the element standard matrices are formed and stored. The INTEGRATION Module embodies the 
integration algorithm for the system of ordinary differential equations as well as access to an algebraic 
equation system solver. For either case, the basic element operation is formation of Eq. (82), and 
assembly of the element matrices into the global representation. The integration algorithm presently 
used by COMOC is an explicit, single-step, multi-stage finite difference procedure with a large region 
of absolute stability (Reference 1 1). A banded Cholesky equation solver is used for the algebraic 
systems. At user-selected points, the OUTPUT Module is called to record the arrays within the 
dependent variable vectors | Q } and other desired parameters. 

NUMERICAL RESULTS 

Numerical evaluation of the finite-element solution algorithm for the two-dimensional Navier- 
Stokes equations, as embodied in the COMOC computer, program, has been performed to assess solu- 
tion accuracy, convergence, and versatility of the finite-element concept. A problem of particular 
value for accuracy studies is developing flow in a rectangular duct of an isothermal, constant-density 
fluid. This problem retains the full nonlinear character of the Navier-Stokes equations, ;and the out- 
flow boundary conditions of vanishing normal gradient for both streamfunction and vorticity. are 
analytically exact for sufficient duct length. Figure 3 illustrates the basic 144 finite element discreti- 
zation of the duct flow problem; for Re = 200, based on duct width, the fully developed velocity flow 
field should be attained within the duct (Reference 13). Additional discretizations studied have 
doubled and halved the number of node rows, and doubled the number of node columns; these are 
referred to, respectively, as the 264, 54 and 276 finite element discretizations. Figure 4 shows the 
COMOC computed velocity distributions at the terminal node column of the solution domain for the 
four studied discretizations. Since velocity is a constant within an element for linear streamfunction 
approximations, Equation (37), the velocity contours are correspondingly constructed. For all dis- 
cretizations, the analytic solution approximately bisects the computed constant velocity within each 
element, and convergence with discretization is illustrated. Within the constraint of element constant 
velocity, solution accuracy is generally good for all discretizations, although refinement of the grid 
certainly improves the solution for points not coincident with an approximate element centroid. 
Poorest improvement in accuracy occurs near the duct centerline where the 2.5% inaccuracy of the 
264 element solution is modestly poorer than either the 144 or 276 element solutions. It is recalled 
that ithe 264 element solution doubles the node rows only, which yields element aspect (length/width) 
ratios of about 50 in this region. Propagation of vorticity to the duct centerline is consequential only 
for x/L > 8.0, and its underprediction in the far downstream region is reflected in the lack of comput- 
ed velocity improvement. This is illustrated further in Figure 5 which presents computed streamfunc- 
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Longitudinal Coordinate - x/L 


Figure 3. Discretization of Rectangular Duct into 144 Triangular Finite Elements 

tion and vorticity in comparison to the analytic values. The improvement in computed vorticity with 
discretization is also illustrated, and the nodal accuracy of the 54 element streamfunction solution is 
remarkable and reflects the well-posedness of streamfunction cast as a boundary value problem. Figure 
6 illustrates computed longitudinal velocity contours at various stations downstream of the duct inlet 
for the 264 element discretization. Velocity overshoot is predicted to occur over the first 4.8% of 
the duct length and is maximum at the 1.5% length station. The magnitude and location of velocity 
overshoot is in agreement with other solutions of the same problem (References 1 2, 14). 

Solution accuracy, as well as computational speed for a transient problem, depends not only 
upon discretization but also upon the manner selected for assembly of the finite-element equivalent 
of the function, g, Equation (82). Initial value integration algorithms are typically derived for or- 
dinary differential equation systems that are explicitly uncoupled in the derivatives; that is, the 
standard form is 

(Q}' = f(|Q|,t) ~ (114) 

Global assembly of the finite element solution algorithm, Equation (82), produces a transient equa- 
tion description wherein the derivatives are explicitly coupled by the [B2,00S] .or [B200SA] finite- 
element matrices. Hence, Boolean assembly of the transient term yields a global description that 
contains a sparse symmetric matrix banded about the main diagonal modifying the { Q } ' matrix. 
The bandwidth is a function of discretization, and premultiplication of the entire equation system by 
its inverse (or equation solver equivalent) is needed to achieve the required form, Eq. (114). The 
inverse is symmetric and a full matrix. Thus, this operation couples the time varying behavior of 
node point dependent variables to all nodes of the discretization. The inversion process can be made 
trivial, and nearest-neighbor transient behavior maintained, by “condensing” the elements of the 
[B200S] and [B200SA] matrices, on rows, which renders the global matrix diagonal. In two-dimen- 
sional space, this operation is a simple averaging, see Table 4. In axisymmetry, the geometry of the 
discretization and lack of translational invariance produces the condensed matrix form listed as 
[B200SCJ in Table 5. 
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Both the consistent and condensed forms for the initial value matrix have been evaluated for 
transient problems by COMOC. Employing the consistent form produces a significantly “stiffer” 
ordinary differential equation system as manifested by the small integration step-size accepted'by 
COMOC. Solution accuracy appears generally acceptable except for spurious behavior during impul- 
sive starts in regions of initially uniform variable distributions. For uniformly positive fluxes into 
these domains, wave-like depression of selective node point values below their initial values is observ- 
ed. Evaluation of this phenomenon, and absolute assessment of solution accuracy was made by 
studying the heating of a quiescent fluid in an a'xisymmetric duct subject to external uniform heating 
from the wall. For this problem, the flow equations are identically satisfied, and only the energy 
equation (70) requires solution after deletion of both the nonlinear source term, Equation (90), and 
the convection term. Figure 7 shows a comparison between COMOC computed radial temperature 
distributions in the fluid, for two discretizations and both assembly options for the initial value 

matrix, and an analytic solution (Reference 15). At the outer wall,-L= 1, solution accuracy for the 

r o 

coarse discretization and consistent matrix assembly is excellent. The depression below the initial 
uniform temperature at the interior radii nodes is illustrated to propagate into the interior region. 

The solution accuracy for the same discretization using the condensed matrix form, is noticeably 
poorer at the wall, and although devoid of the depressions is only a local improvement elsewhere. 
However, doubling the discretization and employing the condensed form produces a comparable or 
uniformly more accurate solution everywhere. 


For small-scale problems involving less than about 50 finite elements, experimentation indi- 
cates that computer CPU execution time for double discretization and condensed form solutions is 
about equal to the corresponding consistent form solution. Even though the consistent form inverse 
can be stored, for larger problems, the lengthy multiplication process coupled with lower stability 
produces execution times uniformly in excess of the corresponding.condensed, double-discretization 
solution. Hence, until integration algorithms for equation systems that are explicitly coupled in the 
derivatives are developed (for example, see Reference 16), it appears uniformly preferable to employ 
finer discretizations and the condensed matrix form for solution of transient problems by finite 
element techniques. 


Abrupt non-uniformities in discretization near a no-slip surface of the domain closure have 
been observed to adversely affected solution accuracy, primarily through the vorticity boundary con- 
dition statement, Equation (103). Evaluation of solution deviations was accomplished using variations 
of the duct flow discretization, Figure 3 after removal of the row of nodes at y/L = 0.95. The analytic 
steady-state vorticity distribution was uniformly applied across the entire duct length as the initial 
condition. The affect on computed streamfunction accuracy, of sequentially replacing nodes along 
the line y/L = 0.95, starting from the left, is illustrated in Figure 8a. At each abrupt change in dis- 
cretization, the streamfunction alternately over- and under - predicts the correct value by less than 
about ± 0.5%. This variation in streamfunction accuracy is probably acceptable for most practical 
flow field predictions. However, near a no-slip wall, these minor errors in streamfunction can 
produce an extremely large error in wall vorticity, as computed through Equation (103) and illus- 
trated in Figure 8b. In all cases, continuing the transient solution from these initial computations 
produced convergence to the correct solution for both streamfunction and vorticity. The direction 
of the diagonals of elements near a wall did not materially affect solution accuracy as long as they 
were all of the same slope, that is, positive or negative. However, Figure 9 shows that alternating 
the slopes has a detrimental effect on computed streamfunction accuracy. As also illustrated, finer 
discretization improved solution accuracy. This evidence indicates that regular and uniform dis- 
cretizations near a no-slip wall are preferable. . 
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Figure 8. Discretization Influence on Vorticity and Streamfunction Near a Wall 
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Streamfunction Deviation 







A crucial analysis for the class of flows governed by the full Navier-Stokes equations is pre- 
diction of regions of recirculation imbedded within an arbitrary flow field. For example, it is 
precisely this phenomenon that prevents use of boundary layer theory in regions of separation in 
otherwise predominantly boundary layer flows. A great deal of laboratory and numerical experi- 
mentation has been conducted for the sample geometry illustrated in Figure 10, which is flow in a 
duct with a considerable enlargement in cross-sectional area. The fluid will not turn the corner, 
but instead provides a smooth area transition by self-generation of a recirculation region. Finite 
difference numerical solution procedures for this problem typically require special handling of the 
corner region to promote generation of the recirculation zone (References 14, 17). There is no 
such requirement for the finite element algorithm and all wall-node vorticities are computed in a 
uniform fashion. Of critical importance for accurate computational representation of this problem 
is prediction of the points of attachment of the streamline dividing the recirculation region from the 
main flow, both at the step and at the downstream reattachment point, points A and B, Figure 10. 
For flow Reynolds numbers based upon step height larger than 20 (Reference 17) the attachment 
at the step occurs just below the corner on the vertical face. The downstream attachment point is 
approximately linear with Reynolds number in the range 20 < Re <200 (Reference 18). The points 
of attachment are identified by a vorticity sign change, and the stream function solution must 
numerically bifurcate to predict recirculation. 
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Figure 10. Flow Over a Rearward Facing Step j 

The finite element solution algorithm, as embodied in COMOC, has been assessed to; accurately 
predict the rearward step steady-state- flow field for Re = 200, treating the problem as transient. 
Shown in Figure 1 1 , is the finite element discretization employed for the problem. The description 
employs 21 1 finite elements and is non-uniform to gain definition in the corner region. Figure 12 
shows computer-generated plots of computed streamfunction and vorticity distributions (upper and 
lower half planes) at various time stations in the transient solution. Figure 12a presents the essen- 
tially inviscid starting solution corresponding to impulsive application of a pressure gradient to a 
quiescent, fluid-filled duct. In Figure 12b, the streamfunction contours have drawn away' from the 
step base and the resultant separation bubble is detected. This bubble rapidly proceeds up the step 
face, Figure 12c, and then proceeds to grow in strength and size, Figure 12d. Figures 12e - 12h pre- 


43 



Figure 1 1 . Discretization of Rearward Facing Step in a Rectangular Duct into 
211 Triangular Finite Elements 


sent growth of the recirculation region to an essential steady-state solution. Note in each Figure the 
predicted occurrenceof a second recirculation region at the base of the step. Within this region, the 
flow field corresponds to flow into a forward facing step, and the fluid again shields itself from 
abrupt changes in cross-sectional area by self-generation of smooth transitions. 

The point of reattachment of the dividing streamline for essential steady state, Figure 12h is., 
about 21 step-heights downstream. This prediction agrees within 10% of an extrapolation of data 
presented in Reference 18 for two-dimensional flows and with the computations of Reference 19 
for axisymmetric flow. No. other, numerical solution to the two-dimensional rearward step prob- 
lem at this Reynolds number is published; further, Dorodnitsyn (Reference 18) observed that 
several independent solution techniques for the finite difference equivalent equations become 
unstable at Reynolds numbers approaching 200. The overall shape of the predicted recircula- 
tion zone, Figure 12h, is in general agreement with the experimental data of Reference 17 ob- 
tained for Re * 1 30. 

Several additional experiments have been performed with the rearward step problem to evalu- 
ate: performance of the versatility of COMOC. Time dependent and non-coordinate surface boundary 
conditions are explicitly acceptable subject to a present limitation in COMOC that rezoning is not an 
automatic feature. Figure 13a shows a close-up of the steady-state rearward step solution, which served 
as the initial condition for the experiments. In the first case, a portion of the duct wall was allowed 
to move outwards in a specified fashion somewhat similar to what might occur for pulsatile flow in an 
elastic vessel. Figure 13b shows the steady-state computed solution after the indentation had been 
stopped and held stationary at 0.4 step heights relative displacement. Computed solution behavior 
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Figure 12a. COMOC Computed Stream function and Vorticity Distribution for Transient Flow 
Over a Rearward Fadng Step. Re = 200, t = 0.000 sec 


Figure 1 2b. COMOC Computed Streamfunction and Vorticity Distribution for Transient Flow, 
Over a Rearward Facing Step, Re = 200, t = 0.200 sec 


Figure 1 2c. COMOC Computed Streamfunction and Vorticity Distribution for Transient Flow 
Over a Rearward Facing Step, Re = 200, t = 0.250 sec 



Figure 1 2d. COMOC, Computed Stream function and Vorticity Distribution for Transient Flow 
Over a.Rearward Facing Step, Re = 200, t = 0.500 sec 



Figure 12e. COMOC Computed Streamfunction and Vorticity Distribution for Transient Flow 
Over a Rearward Facing Step, Re = 200, t = 5.008 sec 



Figure 12f. COMOC Computed Streamfunction and Vorticity Distribution for Transient Flow 
Oyer a Rearward Facing Step, Re = 200, t = 9.998 sec 


Figure 12g. COMOC Computed Streamfunction and Vorticity Distribution for Transient Flow 

Over a Rearward Facing Step, Re = 200, t = 22.365 sec 
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Figure 12h. COMOC Computed Streamfunction and Vorticity Distribution for Transient Flow 


Over a Rearward Facing Step, Re = 200, t = 29.209 sec 




Figure 13a. COMOC Computed Steady-State Streamfunction and Vorticity 
. for Flow over a Rearward Step, Re = 200 


Figure 13b. COMOC Computed Steady-State Streamfunction and Vorticity 
for Flow over a Rearward Step in Irregular Shaped Duct, Re = 200 




Figure 13c. COMOC Computed Steady-State Streamfunction and Vorticity 
for Flow over a Rearward Step with Internal Obstacle, Re = 200 


was stable and the recirculation region has simply readjusted to fill the new interior domain. (The 
indicated appearance of vorticity lines crossing is strictly a function of the high order curve fitting 
operations in the automatic plot package.) The downstream point-of-attachment of the dividing 
streamline was essentially unaffected by the changed wall shape. However, as shown in Figure 13c, 
the insertion of a blunt body into the solution domain drastically alters the predicted recirculation 
region. The inserted obstacle was tapered downstream to preclude the shedding of vortices, hence 
allow attainment of the illustrated steady-state, time-independent solution. The computed vorticity 
solution is markedly altered from Figure 1 3a, and the dividing streamline reattachment has moved 
to within about four step-heights downstream. The second recirculation region at the step base 
is retained in both sample solutions; its apparent unstreamlined shape could be improved with added 
discretization, and resolution could be improved if required. No independent measures are available 
to assess the accuracy of either of these solutions, but the computed solutions do appear acceptable 
on" a physical basis. They amply illustrate the boundary shape variability that is automatically accept- 
able to COMOC. 

Equation (75) has been evaluated for computation of static pressure distributions along the 
centerplane of the duct flow problems. For two-dimensional steady flow with constant density and 
viscosity, Eq. (75) evaluated on the centerplane becomes 

Ap =-A(- \ p , 2 ) 2 /go, 2 dx, +~~ A !/, 12 ^,2 dx, (115) 

p Re J pRe J 

In Eq. (1 15), the first term represents the purely dynamic pressure change due to streamline conver- 
gence. The second accounts for viscous dissipation, and the generic third derivative in streamfunction 
has been replaced in terms of vorticity through Eq. (40). The last term contains the viscous contribu- 
tion stemming from non-parallel flow. Evaluation of Eq. (115) between any two points provides the 
corresponding pressure difference. For fully developed, parallel flow, only the middle term in Eq. 

(1 1 5) is non-vanishing and the expression is analytically exact. For general flows, evaluation of first 
and second order derivatives of tabulated functions is required. In keeping with the use of linear 
finite element approximation functions for computation of the flow field variables, low order dif- 
ference formulas appear appropriate. 

Shown in Figures 14 and 15 are the calculated static pressure distributions along the center- 
plane of the duct and the rearward step problems, for Re = 200, for the computed steady-state stream- 
function and vorticity distributions previously discussed. The total viscous contribution to Eq. (115) 
is separately noted in each case. For both problems, the dynamic contribution dominates the distri- 
bution. For the duct case, the viscous contribution near the centerplane is quite small for about 30% 
of the duct length. For fully developed parallel flow, the pressure decay is solely due to vorticity, and 
the analytic pressure gradient for fully developed flow is also shown in Figure 14. The computed sol- 
ution for viscous contribution to pressure decay and total decay are observed to both be approaching 
this value. Concerning accuracy, the computed streamfunction distribution predicts fully developed 

parallel flow to within 0.75% across the entire duct terminus. The maximum computed deviation 
(~ 25%) from fully developed vorticity occurs adjacent to the centerplane at the duct terminus. The 
error in computed vorticity signficantly affects the accuracy of viscous dominated pressure compu- 
tations which accounts for the slope dispairty at the duct terminus. In most practical flows of in- 
terest, viscous effects are strongest near the containing walls. For the duct flow case, the computed 
terminal vorticity distribution near the wall is accurate to within 2%; a corresponding improvement 
in accuracy of computed pressure in this region would accrue. 
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Solution of the Navier-Stokes equations is typically required for internal- that is, confined- 
flow problems as has been illustrated. A vast number of practical aerospace and aeronautical 
flow field configurations correspond to external flow problems wherein a significant simplification 
to the Navier-Stokes system exists as “boundary layer” theory. This equation system is parabolic for 
steady flows over a wide range of Mach numbers, even for three-dimensional problems, and its solu- 
tion has been greatly explored, see Blottner (Reference 20). However, with the blown flap systems 
on todays aircraft, and/or the supersonic and hypersonic cruise speeds of the space shuttle and hyper- 
sonic research vehicle, there are numerous opportunities to form local regions of imbedded recircu- 
lation within a predominantly boundary-layer-type flow field. Classical examples for supersonic 
flows are shown in Figures 1 6 and 1 7 which illustrate how complex shock wave interaction can in- 
duce local separated flow regions. Such occurrences are of critical design importance since aerody- 
namic control is affected, and local heat transfer rates are significantly increased in these zones. At- 
tempts have been made to use the boundary layer equations for solution of shockwave-boundary 
layer interaction problems (Reference 21), including the neglect of instability-producing negative 
convection terms (Reference 22). These methods have generally proven unsuccessful since the 
numerical mathematics are basically incompatible with the physics of flow. The compression corner 
problem is somewhat more tractable since the approximate location of the shock, hence recirculation 
zone, is at least known. Extension of Navier-Stokes incompressible flow finite difference techniques 
to this problem are reported (Reference 23), and a comprehensive numerical solution of the entire 
transient, viscid-inviscid interaction region flow field is recently published (Reference 24). The rear- 
ward facing cone-step problem is particularly difficult; factors affecting the extent of the recirculation 
region are reported (Reference 25); and the results of a numerical study of the companion problem 
of base flow behind a wedge are published (Reference 26). 

The developed finite element solution algorithm for the steady Navier-Stokes for a compress- 
ible fluid is potentially applicable to this problem class in conjunction with appropriate means for 
establishing the required boundary conditions. For interaction problems, this typically includes 
knowledge of shock impingement location, gs well as inflow distributions and outflow location and 
parameters. Inclusion of density as a computed variable is concomitant with addition of an equation 
of state, Eq. (7). From the definition of stagnation enthalpy, and for a perfect gas, a quadratic ex- 
pression for density is obtained after substitution of Eq. (37) for velocity. Only the positive root is 
admissible which yields the following equation for density in terms of static pressure (p s ), stagnation 
enthalpy (H) and streamfunction derivative. 


P 

The finite element solution requires establishment of density at node points of thb solution domain. 
No complication arises from either pressure or enthalpy, since their nodal values are directly available. 
However, for linear approximation functionals, the streamfunction derivative term becomes 

(*m (*lm [B211S] (117) 

which is strictly element dependent. To establish a unique nodal value of density, Eq. (1 16) is evalu- 
ated within an element loop in COMOC, and nodal values of computed density summed into a global 


c p„ — + 
P Fs R 


c pPs W 
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array. The summation is then divided by the number of entries to establish an averaged value repre- 
sentative of density at a node point. This approach is somewhat similar to certain finite difference 
techniques that employ nodal values for dynamic variables and cell averaged thermodynamic para- 
meters. 
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The utility of the finite element approach for analysis of compressible flows has been evaluated 
in a cursory manner for two problems. Additional downstream discretization was added for the rear- 
ward step configuration, Figure 1 1 , and the solution parameters scaled to simulate isoenergetic flow of 
air al M w = 0.6, based upon a mass-averaged inlet velocity of 122 m/s (400 ft/s), and Re = 200 based 
upon step height. Under these conditions, compressibility and temperature-dependent viscosity (via 
Sutherland’s Law) should be measurable, but their combined effect should not measurably alter the 
essential Reynold’s number invariance of the flow field. As with previous evaluations, the duct flow 
was initiated impulsively and the computed recirculation zone allowed to develop during the transient 
solution. Essential steady-state flow was achieved within 0.15 second; at this time, the downstream 
duct terminus contained fully developed flow to within 0.125% maximum deviation, employing the 
vanishing normal gradient boundary condition for both streamfunction and vorticity. As expected, 
the contours of the computed distributions appeared very similar to those of Figure 12h. The down- 
stream attachment point of the dividing streamline was identically located; the attachment at the step 
had moved to the comer, although a refined discretization would be required to accurately evaluate 
its location for this test case. As before, a second recirculation zone was predicted at the base of the 
step. The null point of the primary recirculation zone occurred on the node column located seven 
step heights downstream of the step. Shown in Figure 18 are computed steady-state distributions of 
the dependent variables across the transverse dimension of the duct at this longitudinal station and at 
the duct terminus. The streamfunction solution illustrates the small flow velocities occurring within 
the recirculation zone; a null point in vorticity is correspondingly indicated. Freestream density 
occurs within this region, and the “pinch” effect of the recirculation zone correspondingly increased 
main flow density by up to 5%. The computed viscosity distributions similarly show largest variations 
outside the recirculation region. 

The finite element solution procedure may be useful for studies of imbedded recirculation 
zones within boundary layer flows, provided boundary conditions are available and shock locations 
known. The supersonic compression corner falls into this category, and a computational solution 
might involve establishing the vorticity-streamfunction equivalent of supersonic boundary layer flow 
as an initial condition, and then bending the plate surface, and skewing the discretization downstream 
of the hinge, to the required angle. Illustrations of sample discretizations for such a procedure are 
shown in Figure 1 9, for the geometry studied by Carter (Reference 24). The existence of the shock 
in Figure 19b could be computationally simulated by the Gump) parallel-velocity boundary condition 
for streamfunction, Eq. (101). As an introductory step to analysis of this problem, the discretiza- 
tion of Figure 19a was employed by COMOC to solve for the vorticity-streamfunction equivalent of 
supersonic laminar boundary layer flow of air for M<» = 3.0, and Re = 7750, based upon distance 
from the plate leading edge to the midpoint of the solution domain. The inflow boundary condition 
for streamfunction, and slope of the top of the solution domain, were established using Eq. (99) and 
a similarity solution for boundary layer flow. The vorticity distribution at inflow was established 
from its definition, Eq. (39), vanished along the top, and was evaluated on the no-slip plate surface 
using Eq. (103). The vanishing normal gradient, although locally incorrect for boundary layer flows, 
was employed downstream for both variables. Shown in Figure 20a is the computed essential steady- 
state vorticity distribution ({ fi)' < 0.001) for the domain of Figure 19a. A high plateau of vorticity 
exists near the plate surface which drops off rapidly as freestream is approached. Figure 20b com- 
pares computed streamfunction and density distributions at x/L = 1 .0 to corresponding values 
obtained from the similarity solution. Agreement is good except for density near the freestream, 
where use of Eq. (116) rounded the distribution. This may result from a poor local speed approxima- 
tion, using Eq. (117), for regions with large gradients in streamfunction; a finer discretization should 
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improve results. The accuracy of computed stream function at the downstream terminus was 
noticeably poorer, as a direct consequence of the applied gradient boundary condition. Judging 
from these results, it appears adequate to extend the solution domain a distance of x/L = 1 down- 
stream of the region of interest for this problem. Furthermore, use of Eq. (116) appears valid for 
compressible flows when combined with attention to discretization. 


CONCLUDING REMARKS 


A finite element algorithm for solution of problems in the two-dimensional Navier-Stokes 
equations has been established. Numerical evaluation within the COMOC computer program system 
has produced a favorable assessment of solution accuracy and convergence. The versatility of the code 
has been illustrated for several problem variations. It appears that the finite element approach of com- 
putational fluid mechanics exhibits much of the required generality and flexibility to serve as an 
effective design and communication tool between the engineer and scientist and the present and next 
generation of digital computer systems. 
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Figure 1 9. Finite Element Discretizations for Supersonic Boundary Flows 
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APPENDIX A 

CARTESIAN TENSORS IN EUCLIDEAN SPACE 


In the rectangular Cartesian coordinate description of Euclidean space, Xj = { x, y, z } , simple 
subscript tensor notation and the summation convention are sufficient to describe the calculus of 
tensor fields. In general curvilinear coordinates, the completely general tensor formalism involving 
co- and contra-variant scalar components and basis vectors is required, see Reference (27). When 
dealing with orthogonal curvilinear (i.e., Cartesian) coordinate systems, simple tensor notation is 
inadequate and the full power of the general tensor is not needed since no distinction exists between 
the two variant forms. 


The appropriate mathematical tool for tensor field calculus operations in Cartesian coordin- 
ates is Cartesian tensor analysis in Euclidean space, Reference (28). The,basic concept of geometry 
is the infinitesimal displacement ds. Denoting yj as the rectangular Cartesian components of an 
arbitrary position vector 7, and xj as the scalar components in any other Cartesian basis, an arbitrary 
vector can be written as 

T= y i ^i = x i e i (A-l ) 

where ^ are the (constant) unit vectors (jj, £) and e- are any other orthogonal basis vector set. The 
infinitesimal displacement (squared) becomes 

ds 2 = d7 - df = dyj&j • dyj a- = dyj dyj (A-2) 


since contraction of Unit vectors in the rectangular basis yields the Kronecker delta.. A similar form 
is obtained for any orthogonal basis, with the exception that, since the basis vectors are not in gen- 
eral unit vectors, their contraction will yield scalar coefficients modifying the Kronecker delta. De- 
noting these as h£, the general form is 


ds 2 = h? dxj dxj 


(A-3) 


where the subscript bar denotes the index is not included in the summation convention. The h^are 
called metric coefficients and are functions of the particular coordinate system-. For example, in 
cylindrical coordinates with Xj = |z, r, 


hj^ = 1 1 » 1 , r [ 


(A-4) 


In spherical coordinates with Xj 


{r, 0, 0} , with <f> the polar angle with domain 0 <x 2 <jr, 



|l, r, rsin0} 


(A-5) 


It is desirable to establish an orthonormal basis for the general Cartesian description. Considering a 
general vector field 7 (as described by the form of equation A-l ), with scalar components vj and 
basis ej, multiply and divide by h: to obtain 




(A-6) 
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In equation A-6, Uj (= hj Vj) are called the physical components of the vector v", and ej are unit vectors 
of the orthonormal basis. 


The fundamental calculus operation requiring identification is differentiation. Referring to 
equation A-6, the derivative of v" with respect to the Xj variables is, using the chain rule 


bv _ b__ 

d X: d X: 


(Ui^j) 


3u: be-. 


(A-7) 


The last term in equation A-7 describes the change in the unit vector with respect to Xj. Since 
this is also a vector, it can be described by scalar components projected on the basis. The coeffi- 
cients of proportionality for this description are the physical components of the Christoffel symbols 
of general tensors, contracted over index pairs. They have been termed “wryness” coefficients, Ref- 
erence (29); signifying them by < 0 ^ and expressing v’ in the basis, equation A-7 becomes 

dv” A A 

- = u i;j^i = (u i,j + cu ijk u k)ei (A-8) 

In equation A-8, the spatial (vector) derivative of a vector is signified by the semicolon. The partial 
derivative of a scalar (component) is denoted by the comma. Since equation A-8 is a vector identity 
the scalar components are identical. Hence 


1 ;j 


- u 


i,J 


+ w ijk u k 


(A-9) 


defines the vector derivative of the scalar components of a vector (first order tensor) in an arbitrary 
orthonormal basis. Equation A-9 can be readily generalized to the vector derivative of the scalar 
components of an n** 1 order tensor, see Reference (28). 

The properties of the wryness coefficients are readily determined. They are not scalar com- 
ponents of a third order tensor since they do not obey the tensor transformation law. All components 
vanish for the indices distinct. The coefficients are skew-symmetric in all index. pairs. Hence, the 
only nonvanishing terms are of the form 


w kjk w jkk 

In particular, the singular nonvanishing coefficient in cylindrical coordinates is 
_ 1 

<*>2 3 3 ~ ~ 

In spherical coordinates, obtain 

_ 1 

0)122 ~ ^'133 jr 


(A-10) 


(A-ll) 


(A- 12) 
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From equation A-9, the vector derivative of a scalar is simply the partial derivative with 
respect to Xj, commonly referred to as the gradient. From equations A-l , A-3, and A -6, obtain 


90 _ A A _ 1 00 A 


(A-14) 


Equation A-14 defines the partial derivative, signified by the comma. Equating scalar coefficients, 
obtain by components 


< >-i “irjb; ( > 


(A-l 5) 


The operation of curl is the contraction of equation A-9 with the general alternating tensor ejj^, 
which is skew symmetric in all index pairs and has positive scalar components for (i, j, k) in the 
cyclic order (1, 2, 3). Hence, 

curl v’-£j = e ijk u k ; j (A-l 6) 

The general alternating tensor is related to the Cartesian alternating tensor tjjk, with unit magnitude 
scalar components, by the determinant of the space metric, which for Cartesian tensors is simply the 
product 

J 4 - hi h 2 h 3 (A-l 7) 

Equating scalar coefficients in equation A-l 6, and using equation A-17, obtain for the i** 1 component 
the curl of the vector "v 


curlv.^ = }e ijk u k;j (A-l 8) 

Equations A-3, A-9, A-10, A-l 5, and A-l 8 are sufficient to facilitate the calculus of Cartesian 
tensors in Euclidean space. Integral theorems, developed for rectangular Cartesian tensor fields, 
transform directly to general Cartesian tensor fields. Of particular interest, the divergence theorem 
becomes, for an arbitrary order Cartesian tensor field T^. . . , 

/ R CTjjk •••>;{ dr - <f dR <Tjj k . . . ) • » 8 d a (A-l 9) 

where is the unit outward pointing normal vector to the closure dR of the domain R. Similarly, 
for Stokes’ theorem obtain 


/ (curl T ijk • • • > * da = ^ dL ( T ijk/ ••)•<£ ds 


(A-20) 


where 3L is the bounding curve of the closure 0R and tj is the unit vector everywhere tangent to 3L. 
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APPENDIX B 

OBSERVATIONS ON THE FINITE ELEMENT SOLUTION ALGORITHM FOR THE 
NAVIER-STOKES EQUATIONS USING LINEAR NATURAL COORDINATE 
APPROXIMATION FUNCTIONALS 


The Navier-Stokes equations, governing the two-dimensional flow of a compressible viscous 
fluid in rectangular Cartesian space, contain several nonlinear terms involving both density and 
viscosity gradients. The finite element solution algorithm retains the vector differential operations 
within its formal definitions, and options exist to perform integrations by parts as well as employ 
integral calculus theorems to simplify terms. This appendix presents observations made during the 
theoretical exploration of the consequence of such operations and their impact upon the numerical 
implementation of the finite element theory. 

General Observation: The finite element solution algorithm applicable to elliptic partial differential 
equations is of the form: 


/ w ^kk dr = / 


Wf, (0*, 0? k > x i> dT + wf 2 (0*) da 


3R 


m 


(B.l) 


where 0* is the approximation to the dependent variable, fj are specified functions of their arguments 
and contain no second derivatives of 0*, and the weighting functions W are identical to the approxima- 
tion function for 0*. Within a finite element, 


^m = |M T M m 

W = | L | 

and in two-dimensional rectangular Cartesian space 




-1 ] 


x 3 - x 2 

i i r 

A 1 


L >k = - 

1 

A * 

• e, + 

- x 3 

l J K X 2 

1.0. ] 

x 2 y 3 

+ x 2 


where 6: are unit vectors in the local reference system, Figure B-l . 


(B.2) 

(B.3) 


(B.4) 


(B.5) 
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Figure B-l. Two-Dimensional Finite Element 

Of particular interest is determining the behavior of Eq. (B. 1 ) when integration by parts is 
performed on terms off, as well as0 tkk • As an introduction, integrate the left side of Eq. (B.l) by 
parts to obtain: 




W 0 , k n k d o-J W, k <£, k d it 

Rr 


‘m 


(B6) 


Theorem 1 : If 0 is a linear function of its arguments, Eq. (B.6) vanishes identically for evaluation 
within a single finite element, 

Proof : 


1 dr = 0 


/ w **k dr = / 1 L 1 { L } T .kk dT =/ 1 

^m ' 


I 1 / 1 

- 1 1 

1 

x 3 -x 2 

L 1 — 

i 

+ 

-x 3 

1 j \x 2 

0 1 

. x 2 y 3 
’1 

X 2 


(B.7) 


£ W (l) , k n k d a = J (l) j l} * nj^da 1 + J {l} {l} T k n" do + f {l} |l) , k n” 1 da 

9R m dRj 3R n 3R in 
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{ ((6 |j|h. , .°i 
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(x 2 y 3 ) 2 


x 2 y 3 


x 3 -x 2 

-x 3 
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.{x 3 -x 2 , -x 3 ,x 2 } ) dr H 
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[B21 IS] m j^) 


m 


(B.9) 


Combining Eq. (B.8) and Eq. (B.9), the identity in zero is accomplished, Q.E.D. 

Theorem 2: If <p is a linear function of its arguments, the surface integral contributions resulting 
from integration by parts, Eq. (B.l), can be made to cancel in pairs upon assembly, by Boolean 
algebra, of the global solution algorithm. 

Proof: Along the line common to the finite elements defined by domains Rf and R jj , Figure B-2, 
specify the local value of the normal gradient of the function <t> as an algebraic average. The coin- 
cident portion of the closed surface integral for each domain becomes 
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Figure B-2. Adjacent Finite Elements 
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Upon Boolean assembly of Eq. (B. 1 2) and (B. 1 3), the first of each equation' is loaded into the global 
vector on $j at the location <j> j. Hence, cancellatipn by pairs is achieved for the surface integral 
contributions stemming from Eq. (B.6), Q.E.D.. 

Theorem 3: If0 is a linear function of its arguments, and Eq. (B.l) involves a scalar nonlinearity 
that is not differentiated, the surface integral contribution stemming from integration by parts can 
be made to vanish in pairs upon Boolean assembly of the global solution vector. 

Proof: The compatibility equation for the two-dimensional Navier-Stokes equation can be integrated 
by parts, as 

/ **>* * = 4 V/j *. k n k ta- f W, k i*, k dr (B. 14) 

R dR P R P 

Using the averaging concept developed in the previous proof, obtain from Eq. (B.l 4) on an element 
basis: 



( ) da = 



(B.l 5) 


Along the line connecting nodes 1 and 2, Figure B-2, Eq. (B.15) becomes, using Eq. (B.4) and (B.5) 
within elements Rj and Rjj, 


f ( ) da 

9Rj 


f ( ) da 
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(B. 17) 
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The only elements of the finite element matrix (—1 involved in non-zero products in Eq. (B.l 6) 

'P *m 

and (B.l 7) are for node locations on the common line. Hence, upon Boolean assembly, the con- 
tributions from surface integrals involving a scalar nondifferentiated nonlinearity cancel by pairs, 
Q.E.D. 


Observation: The convection term, involving the curl of the streamfunction (vector potential 
function), appears in the differential equation and finite element description for each dependent 
variable except streamfunction. An integration by parts is possible. 


Theorem 4: Within the convection operator, and for <p a linear function of its arguments, no 
arbitrariness exists for evaluation of the generated surface integral. 
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Proof: The typical convection term is, upon integration by parts, 


/ We 3ki (co<//, i ), k dT = j> We 3ki co\[/, i n k da- J W, k e 3ki wtf.j dr (B.18) 


R 


m 


3R 


m 


R 


m 


For the surface integral, using Eq. (B.2) — (B .5), obtain 

/ We 3ki coV/,j n k do = f |l| {l} T { ft| m e 3ki (4% n k do M m (B.19) 


dRm dR m 


Note that the gradient operator on | L } is orthogonal to the local normal via the alternator. From 
the natural coordinate function description, Eq. (B.4), such derivatives involve end point nodal values 
only. Hence, for integration over the three sides of the finite element, Figure B-l , obtain 
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(B.20) 


Evaluating the remaining term on the right side, and subtracting it from Eq, (B.20) produces the 
identity. Hence, there is no arbitrariness in evaluating the terms stemming from this integration by 
parts, Q.E.D. 

Straightforward evaluation of Eq. (B. 18) for linear functionals yields 
/ We 3ki ( W *, j ). k dr = / |L). 3kl MllL}, ie3ki |L) T k dxln} (B.21) 
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The alternator contraction produces the skew-symmetric matrix. 
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Equation (B.21) takes the form: 


/ We 3ki dr 



-1 

0 

1 



(B.22) 


(B.23) 


Observation: The finite element solution algorithm produces several nonlinear matrix forms that 
are eligible for integration by parts. For example, in the differential equation for vorticity, an 
additional convection term results for nonuniform density flows of the form: 


67 



J ^ c 3ki 2 

This term vanishes identically, for linear functional approximations, unless an integration by parts 
is performed. Hence, 

/ ( )dr= ^ <“>•/ w .k'3ki(;).k{«'.S> !d ' 


(B.24) 


Equation (B.24) vanishes identically for \p,% evaluated solely within R m . However, 

Theorem 5: Using linear functional representations for all dependent variables, the surface integral 
contribution to Eq. (B.24) can be made to cancel by pairs, upon Boolean assembly of the global 
solution algorithm, by averaging the available derivative operations over closure segments of the 
respective finite elements. 

Proof : Referring to Figure R-2, along the common closure segment connecting nodes 1 and 2, the 
surface integral contributions within Rj and Rjj can be written as 

x 2 

/ ( )da = i / (Lj W?fL I },(L n i,lW II e 3ki {LL k T (I}n! dx (B.25) 

dRj o 1 

/ ( )da= T / MMii fLlll.el L ll.g MlM* Clx (B - 26) 

3R n o 

Since the inner product on i //,£ is commutative, and since the alternator e 3 ki removes arbitrariness 
in the evaluation of |l| ^ , and since nj are antiparallel, then upon Boolean assembly of Eq. (B.25) 
and (B.26) into the global vector, cancellation by pairs is achieved, Q.E.D. 

Completion of the matrix form of Eq. (B.24) is straightforward. 

/ ()*•■-/ W, ie 3 ki (i). k i(*, s )=dr 

Rm Rm 


The inner contraction upon i//,g involves the matrix [B21 IS] , Eq. (B.8), while the alternator con- 
traction produces the negative of Eq. (B.22). Hence, the final form becomes: 


/( )dr 
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+ ± {*} T [B211S] { 4 ,} 
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(B.27) 
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APPENDIX C 


SYSTEM INVARIANCE UNDER COORDINATE TRANSFORMATION 

Use of a locally defined coordinate system for matrix moment evaluation is desirable from the 
accuracy standpoint and mandatory for evaluation of the nonlinear terms in the Navier-Stokes equa- 
tions. The classical Method of Weighted Residuals (MWR) solution to a differential equation, N(q) = 
0, involves approximation of the dependent variable by a function, q*, and rendering the error ortho- 
gonal to the solution by formation of the residuals 



The finite element algorithm transforms formation to a finite element R m , and there are as many equa- 
tions (C. 1 ) as unknown expansion coefficients in q^, on a local basis. The theoretical complication 
stemming from gradient boundary conditions can be neglected for this development. 

Use of natural coordinate functions in a local Cartesian basis is preferred; hence, it is necessary 
to transform Eq. (C.l ) from the global basis of its derivation. Assuming use of truncated power series, 
and to within an arbitrary coordinate translation, see Figure C. 1 , the global reference approximation 
function can be expressed as 

q*( Xi ) = tx-X]T[r] m (C.2) 

In Eq. (C.2), the elements of the coefficient matrix [r] m involve only node point coordinates of R m 
and are readily evaluated (reference 12). It is required to establish Eq. (C.2) in the x- basis, related to 
Xj by the tensor transformation law 

xj = 0y (xj-Xj) (C.3) 

or equivalently, in matrix notation 

|x'} = [0] jx-X} (C.4) 

In the primed coordinate system, Eq. (C.2) takes the form 

s l x ’l T t r 'lm iQlm <C.S) 

Combining Eqs. (C.4) and (C.2), and comparing to Eq. (C.5) provides the definition of [T'] m as 

ir'] m = [0 ] " 1 [r] m (C.6) 

Identifying the weighting functions for the Galerkin criterion within MWR as 

w m == [r] T Ix-x} (c.7) 

J m 
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= x , sinx 2 dr 

(C.l 1) 

For cylindrical coordinates and' axisymmetry, 


x 2 = + x, sin a + x 2 cos a 

(C.l 2) 

where x 2 ^ ) is the radial coordinate of node 1 for R m , Figure C.l . 
natural coordinate functions defined by 

The equivalent expression for 

<lm = )L(x;)| T 

(C.l 3) 

becomes 


x 2 = x 2 ^ ) + x- Lj 

(C.14) 
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